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1 Introduction 


1.1 Background 

Since the early days of artificial satellites it has been recognized that one could map varia- 
tions in the Earth’s gravitational field by accurately tracking the perturbations in satellite 
trajectories. Unfortunately, however, the accuracy limitations of conventional absolute 
tracking methods make these single satellite methods inadequate for high resolution map- 
ping. A way out of the impasse is to use the so-called satellite-to-satellite (SST) tracking 
methods. The first of these is the high-low system in which a test satellite orbits in low 
earth orbit where its trajectory is most affected by the disturbing geopotential and use one 
or more satellites in well known orbits at higher altitudes to track it. The method which 
we are studying here is the low-low, interferometric, or differential method. In the ideal 
configuration, a pair of satellites orbit in identical low polar circular orbits, one trailing the 
other at a “fixed” distance. Millimeter or optical electromagnetic waves are transmitted 
from one to the other and back so that the relative speed of the two satellites can be 
detected as a doppler shift in the wave frequency. Our area of interest lies in developing 
the methods by which one processes this relative range rate data to map the geopotential 
field. 

Important work in this area has been done by Colombo [4], who proposed a method 
based on random fields on a sphere. See also our own work, Bose et al. [1], in which 
the advantages of treating the global geopotential as a homogeneous and isotropic random 
field on the sphere and sampling it on a uniform global grid are analyzed in detail. In a 
later work, Colombo [6] presented the theory and results of a variational method motivated 
by G. W. Hill’s lunar theory. One of the conditions for Colombo’s methods is that the 
satellite pair execute synchronous orbits, i.e. that the number of sidereal days and the 
number of revolutions be relatively prime causing the trajectory to close on itself and the 
ground track to repeat. 

The reader is also referred to Wagner [16], who presents spectral analyses based on the 
simplistic signal equation of Wolff [17], and performs a low degree and order 4x4 recovery 
based on a more elaborate signal equation derived from conventional Keplerian element 
perturbations. 

Most important as a basis for the present work is Kaula’s paper [10] and the subsequent 
work of Bose and Thobe [2] in which an observation equation in the frequency domain 
was developed starting with a harmonic expansion of the geopotential in terms of modified 
Keplerian elements of the satellites. Fourier analysis of the harmonic expansion leads to an 
approximate uncoupling of the signal equation according to the harmonic order and parity 
of degree. Included in [10] are Kaula’s results of numerical recovery experiments for an 
8x8 degree and order harmonic geopotential. In Thobe et Bose [14], we applied Colombo’s 
condition of a synchronous orbit to the method of Kaula [10], rigorously demonstrating 
the uncoupling of the signal equation and the resulting mathematical simplification. 
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1.2 Purpose and Scope 

The purpose of this research is to develop techniques of recovering the spherical harmonic 
coefficients of the geopotential field from SST range rate data, and to demonstrate their 
effectiveness by numerical experiments. A key objective is to improve the resolution of 
gravity models, which implies recovering harmonics of high degree and order. At the 
same time, adequate numerical accuracy and reasonable computational efficiency must be 
maintained. Care had to be taken at all stages of development to insure that our soft- 
ware would be able to deal with much larger models than most previous efforts, nominally 
180 x 180 and potentially still larger. The analysis had to be rigorously correct yet simple 
and elegant so that exploitable regularities would be apparent, and the equations could be 
reduced to their simplest form. Data storage requirements had to be minimized. Efficient 
methods of computing special functions had to be devised. Observation matrix struc- 
ture had to be exploited to reduce computational complexity to a manageable minumum. 
Numerically stable solution techniques had to be selected so that mathematically correct 
solutions would not be swamped by numerical errors. 

1.3 Technical Approach 

The first step toward achieving the objectives set out above was to rederive the signal 
and observation equations from fundamental physical principles. The rotating earth gives 
rise to a time dependent gravitational potential, which in turn means that the mechanical 
system of a free particle, i.e. a satellite uninfluenced by drag or other external forces, can 
be characterized by a time-dependent hamiltonian. By using the basic variational calculus 
energy conservation laws as they apply to a time varying hamiltonian, and applying the 
elementary relationship between speed and kinetic energy, we were able to derive a simple 
yet correct signal and observation equation. We thereby avoided the complicated construc- 
tions from perturbations of individual Keplerian elements characteristic of previous signal 
equation developments. 

The next step was to impose the condition of a tuned synchronous orbit of N T revo- 
lutions in Ni days where these are relatively prime integers. Geometrically, this means 
that the spacecraft covers the earth’s surface in a set of evenly spaced north-south and 
south-north traverses and repeats the same trajectory relative to an earth-fixed frame once 
in each mission period. Mathematically, it means we are dealing with a periodic process 
and so can apply old fashioned Fourier series. 

After applying the Fourier coefficient operator to the time domain signal equation to 
obtain the SST range rate spectrum, we find that all observation matrix elements vanish 
except at the zeros of an integer relation connecting the frequency k, and the harmonic 
indices /, m, p, and q. This gives rise to a sparse structure — in fact, it uncouples the 
signal equation into smaller independent sub-equations. With such a sparse structure, the 
previously massive computational requirements are reduced by several orders of magnitude 
to become manageable on a modest computer, such as the DEC VAX 11/750. By factoring 
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the observation matrix into the product of a diagonal matrix factor, a real matrix factor, 
and another diagonal matrix factor, and employing an inclination function identity, we 
reduced the observation matrix to a real one with half the number of rows as the for- 
mer complex observation matrix. Remarkably, the reduced matrix is simply an array of 
inclination functions. 

As to software, we set it up to perform the Fourier analysis on the entire mission 
range rate history at one time using a public domain mass storage Fast Fourier Transform 
(FFT) package adapted by us to handle double precision computer arithmetic. We then 
use sorting to group elements of the spectrum into their respective blocks. Thereafter 
each block is processed separately: the matrix elements are computed and the observation 
equation is solved. In order to demonstrate the essentials of the method, we made certain 
interim simplifications: we took the reference orbit to be circular, neglecting linear and 
higher order terms in the eccentricity. Also we limited ourselves to a single iteration in the 
solution of the observation equations. 

Key to being able to compute the range rate spectrum and matrix elements, was the 
invention of an efficient, numerically stable, partially recursive technique for computing the 
inclination functions. Naive application of published analytical formulas led to problems 
of poor efficiency or instability. 

The experimental phase of the study began with the development of software to perform 
the Fourier analysis, to compute the inclination functions, to construct the observation vec- 
tor and observation matrix block by block, and to solve the observation equation. Careful 
testing and validation of each element of the software was a necessity. An analytical sim- 
ulation totally under our control was developed for debugging and preliminary validation 
of the software. Finally, we tested the method and software on foreign simulation data 
from the University of Texas Encke method numerical integration of a 32 sidereal day, 525 
revolution synchronous orbit mission. 

1.4 Summary of Results 

We analytically simulated a one-day, 16-orbit mission, producing artificial range rate data 
from which we successfully recovered the geopotential coefficients to an accuracy of ap- 
proximately one part in 10 4 . We then tested our method and software on data produced 
at the University of Texas by Schutz et al. [13]. Using the SST range rate signal from 
their much larger scale 32-day 525-revolution simulated mission, we computed the com- 
plex fourier coefficients of the range rate signal, evaluated the observation matrix elements, 
and numerically solved the observation equation for several of the lower degree blocks to 
order 180. For instance, for order m = 1 and even degrees l < 40, errors in the recovered 
coefficients were typically as low as 10%. The error figures deteriorated rapidly for higher 
orders, unfortunately. We ascribe this in part to fact that we assumed a very simple 
circular reference orbit, neglected linear and higher terms in the geopotential harmonic 
expansion, and did not iterate our solution. On the other hand, the spectra predicted by 
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our signal equation show good qualitative agreement with the spectra we extracted from 
the Texas data. This can be seen by visually comparing the spectral magnitude plots in 
the Appendix A. It is possible that errors in the Texas simulation contributed in part 
to the discrepancies; in fact, the spectra of the Texas data show unnatural regularities at 
wave numbers greater than 50000 (half the presumed bandwidth), tending toward constant 
amplitude and phase. 

A gratifying conclusion to be drawn from these experiments is that high resolution 
recovery can be performed on modest sized computers. Each block takes approximately 5 
minutes on a DEC VAX 11/750. Further work is needed however to improve the quality 
of the recovered coefficients. 

1.5 Overview of Report 

In this report we present both the theory and validating numerical experiments of our 
own system of geopotential coefficient recovery, as it has evolved from the seminal work of 
Kaula. In Section 2 we derive our new signal model starting from fundamental physical 
principles. We first obtain a spherical harmonic expansion of the kinetic energy of a single 
satellite, then its speed, and then the relative speed of two satellites. This is in constrast 
to previous methods, which either simplistically relate the relative range rate to potential 
energy differences between the two satellites, or employ complicated constructions based on 
perturbations in individual Keplerian elements. In Section 3 we develop the signal model 
into an observation equation which can be solved to recover the coefficients. The sparse 
structure of the frequency domain observation matrix is exploited to drastically reduce the 
cost of obtaining a solution. We take advantage of additional, newly discovered symmetries 
to reduce the observation matrix still further, finally obtaining a block diagonal matrix 
of real trapezoidal blocks. In Section 4 we present certain numerical techniques crucial 
to generating the matrix elements and solving the observation equation. In particular, 
we present our partially recursive, summation free method of computing the inclination 
functions Fi mp (i). Also we outline the use of numerically stable Givens orthogonal rotations 
to obtain the optimal solution of the over-determined linear observation equation. In 
Section 5 we describe our validation methods, especially the analytical simulation which 
we used to produce test data. Then we detail the results of large scale (maximum degree 
and order 180) coefficient recovery experiments from the University of Texas 32-day, 525- 
revolution Cowell simulation. Section 6 presents our conclusions and recommendations. 
In Appendix A we have reproduced tables and spectral plots showing partial quantitative 
results of our recovery experiments from the Texas data. 
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2 Signal Model 

2.1 Single Satellite Motion 

The geopotential is given in Kaula [9,10] as the spherical harmonic expansion 

^-r + EE'i. (1) 

' 1=2 m=0 

where in terms of the geocentric latitude <j> and longitude A, an individual harmonic may 
be expressed as 

tm l 1 00 

Vim = - sin (i>){Clm COS m\ + Sim sin m\) = Y,Y, V ‘mpq ( 2 ) 

' p=0 <j=— oo 

In terms of the Keplerian elements of an orbiting satellite, and the sidereal time 6 expressed 
in radians, the potential is expressed as 

Vimpq = —^j^Fi mp (i')Gip q (e)ReSimp q (u;,M,Q,,8) (3) 

Note that the physics sign convention is used for the potential energy V, whereby a = 
— VV. Here a denotes specific force or acceleration. Other symbols in Equation 3 are 
given by 

Slmpq ~ -dim GXp ilfilmpq (4) 

i’lmpq = (l ~ 2 p)u> + (l -2p+ q)M + m(Q - 6 ) (5) 

Alm=(-i) ll - mU [Clm-iSlm] (6) 

The notation \l — m| 2 is equivalent to l — m mod 2. The cosine and sine harmonic co- 
efficients Cim and Sim, respectively, are a set of empirically determined numbers which 
quantitatively describe the disturbing potential. The objective of this study is to develop 
effective techniques for estimating these coefficients from experimental data, namely the 
satellite-to-satellite (SST) range rate A R. In the following, we shall estimate the complex 
coefficients Ai m and then solve Equation 6 to obtain 

Cim = ReiV-”\>Ai m (7) 

Sim = ImtH-"*M /TO (8) 

In addition, the special functions of orbital mechanics, the inclination function Fi mp and 
the eccentricity function Gi pq , have more complicated definitions for which the reader is 
referred to Kaula [9]. Techniques for computing these functions are discussed in Section 4.1. 

We introduce the total energy E, the kinetic energy T, the potential energy V , the 
speed or path rate i, the partial time derivative d/dt, and the total time derivative d/dt 
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or overdot (’). The hamiltonian H is an explicit function of the position x 

\ the momentum 

p, and the time t: 

E = T + V =\s 2 + V = H(x, p, t) 

(9) 

From the fundamental variational theory of classical mechanics we have: 

• dH 
E= 9t 

(10) 

From Equation 9, we have 

dH _ dV 
dt dt 

(11) 

By definition, 

dV 
V — — 

dt 

(12) 

so that 

■ p ^ dV dV 

T ~ E V ~ dt dt 

(13) 


Integrating Equation 13 and applying the perturbation operator A to Equation 9 we get 
a general expression for the perturbed speed of a satellite. 


As = 


f [dAV 

dAV 

/ [ dt 

dt 


dt 


(14) 


We observe that of all the dependent variables in Equation 3, only 9 is explicitly a function 
of time. The others collectively represent the position x and momentum p of the restricted 
phase space. We need not be working in a canonical coordinate system for our derivation 
to be valid; we only require that the time dimension be distinct from the position and 
momentum. To compute T using Equations 3 and 13, we need the time derivatives 

9 R gf O', M, a, u».lt - to] + »o) (15) 


d Re Si 


mpq 


dt 


(uj, n[t — to] 4- Mo, fi, w e [i — to] + #o) 


(16) 


There is a degree of approximation in Equation 16 because uj, n, M 0 , and Q have been 
treated as constants. A better approximation would incorporate at least the linear terms 
in t due to secular perturbations. Continuing: 


dSlmpq _ . .dlpi 
— A i m Z 


dt 


mpq __ 


dt 


— iftlUJ^Slmpq 


= = •’[('- 2p +«)» - "■".I Sr 


mpq 


(17) 

(18) 
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d R.C Slrnpq 

dt 


d Re S l mpq 

dt 


— Re l{J> 2p “f" q^TlS 'irnpq 

— (I 2 p (^)n Im Siftipq 

Substituting into Equation 5 the relations 

M — Mo -f- M[t — to] 

0 — 9q + 0[t — to] 

M = Nri> 

8 = u e = N d ip 

ifrimpqi 0) = {l — 2 p)u + (l — 2p + q)Mo + m(Q — 8o) 

V’(t) = i>[t - to] 

we get an expression for which is linear in time: 

</> lmpq(t ) = *l>lmpq(0) + [(/ “2 P + q)N r - mN d ]lp(t) 


Finally, we get the kinetic energy perturbation 

A Ti mpq = J 


ATlmpqdt = ^F lmp G lpq 


a 


f+i 


[l-2p + q]N r Re Si mpq 

[/ — 2p + g]JV r - mNd 


(19) 

( 20 ) 


( 21 ) 

( 22 ) 

(23) 

(24) 

(25) 

(26) 


(27) 


(28) 


The kinetic energy perturbation is related to the velocity perturbation as follows: 

AT = A (|u -vj = v- Av (29) 

which is to say that AT is proportional to the along track component of the perturbation 
in velocity. Applying ««na for a circular orbit, we get: 


_ /j,a l e ~ [l — 2p + q]MKeSi mpq 

[;_ 2 p + q ]M-m» 


(30) 


Thus we have expressions for the perturbed kinetic energy and speed of a single satellite 
in a circular orbit. 

We note that for a circular orbit the angle between the satellite-to-satellite relative range 
vector R and the velocity vector v is \8u, which is rather a small angle. For computing the 
range rate perturbation, we are interested in the component of Av on R. It is not quite 
accurate to substitute the one component of Av for the other, but this is what we have 
done, and the results of this approximation are reflected in the observation Equation 43 
and in the numerical experiments reported in this document. A more rigorous observation 
equation derivation, which does not make this approximation is carried out in Section 2.3. 
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2.2 Relative Motion 


In this section we apply the results of the previous section to compute the perturbed 
relative speed, a.k.a. the perturbed SST range rate, of the two satellites. We introduce 
the difference operator S(-) = (•)! — (-) 2 , i.e. the value of some quantity for the number 
1 (lead) satellite minus that for the number 2 (trailing) satellite. Since the orbits are 
coplanar we have: 


8St mpq = Ai m exp *[(/ - 2 p)ll> + m(fi 


x {exp i(l — 2p + q)(M + S -^~) - exp i(l - 2p + q)(M - 
= 2 ■ sin (1 -2 p + q)^- 


(31) 


(32) 


and 

£ 

8 Re Si mpq = -2 sin(Z -2 p + q) — • Im Si mpq (33) 

From the geometry of the circular orbit we obtain the SST range rate perturbation in the 
time domain: 


. • ... oivi , . 

^ Rlmpq = uAsimpq COS “ (34) 

n 8M . JM T „ , , 

= -2 cos — sin(/ -2 p + q) — Im Si mpq x (35) 

r ~ 2 1± ^ t 

na i +2 ,m P lpq [I - 2p + q]M - m6 

Here M = (M\ + M 2 )/2 is the mean anomaly of a fictitious mean satellite. Next, we 
compute the fourier coefficients of the SST range rate Equation 35. Applying Equations 4 
and 27 we get 

SimM) = Si mpq ( 0) exp i [(/ -2 p + g)N r - mN d ) ip (36) 

=v =w 

Applying the fourier coefficient operator T , defined as 



F[f(4 0 1 k ] = ^ J Q exp (-ikxj))f(4>)d4> 

(37) 

we have in general 




T[ Im(u exp iw\l>) | fc] = — {8 w _ k v - 8 w+k v} 

Zi 

(38) 

where 8 n = 1 if n = 

0 and 8 n = 0 otherwise. In particular, 


f[Im Si mpq (iJ>) | 

^ ” ^{$(l-2p+q)N I -mN <i -kSlmpq(0) ^(f— 2p-\-q) N T — mNa+k *$7mpg(0) } 

(39) 
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from which we obtain the frequency domain expression for the SST range rate: 




pal 01V1 ... . orn [i -zp-\- q\l\ T . . 

cos T Sm( “ P + q) ~2T ' [i-2p + g)W r -mJV d X(40) 

^ { < ^(Z-2p+g)N r -m./V d -A:‘S'/TOpg(0) ^(/-2p+g)Nr— m.N d +fc , S/mpg(0)} 




/M [/ — 2p + q\Nr 


For the sake of simplicity, let us consider only the case m ^ 0. Furthermore, let us 
consider / < / max where / max is a suitably chosen function of N r and jVj. Then, because of 
redundancy, we need only consider the possibility 


k = (/ — 2p + g)iV r — mNd (41) 

This allows us to drop the terms corresponding to 

— k = (/ — 2p + g)JV r — mNd (42) 


to produce the simpler SST range rate equation: 


^<4 


^Rlm’pqik^) — ® COs(JVi 


na 


sin (i l ~ 


, 6tl>' 

2 P + q}N r f 


[/-2p + g]iy r 5 , / mpg (0) 

[/ - 2p + 4 ]JV r — mNi 

(43) 


We arbitrarily exclude the case k = 0 because of the vanishing denominator. We have 
reached our goal of a signal equation in the fourier domain. 
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2.3 Extended Signal Model 

In the previous section, we derived an SST range rate perturbation signal model based 
on the kinetic energy perturbation of each of a pair of satellites. The kinetic energy 
perturbation is easily converted to a speed perturbation which is essentially the horizontal 
component of the velocity perturbation. This is projected onto the satellite-to-satellite 
range vector. What is neglected in this procedure is the vertical component of individual 
satellite’s velocity perturbation which also has a non-zero projection on the range vector. If 
we pass to the limit as the satellite separation distance approaches zero, clearly projection 
of this vertical component on the range vector approaches zero as well. However, numerical 
experiments reported by Colombo [], indicate that for reasonable separations, the vertical 
components can contribute significantly to the overall SST range rate spectrum. 

In this section we present an extended range rate signal model which includes the 
vertical components. The perturbation of the vector velocity of a satellite cannot be 
obtained from consideration of a scalar kinetic energy, so the approach will be different. 
However, the result is essentially in agreement with what we have already derived except 
for the addition of a term proportional to the sine of half the angle separating the two 
satellites. 

Following Bose and Thobe [2, sec. 4] we compute the inter-satellite range 


R=\\R\\ 

(44) 

the range rate 


R — R • R/ R 

(45) 

and the range rate perturbation 


AR = R- 1 {£ • A^ + R ■ [l - R~ 2 R ® £] • A 

(46) 


Neglecting the term containing the projection operator 1 — R 2 R® R and converting to 
the polar coordinates r and u in the orbit plane, we get 


. • . OU ... . . . OU ... . . . , 

AR « sin — (Ari + A r 2 ) + r cos — (Aui — Au. 2 ) (47) 

Zd A 

or in terms of the operators S and 8 defined by 

*(•) = (Ox - (■)> (48) 

and 

SO = (•)!+(•)> (49) 

where 1 and 2 represent the leading and trailing satellites, respectively, we get 

AT? « sin ^E( Ar) + r cos ^8(Au) (50) 

£ Z 
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We again use the physics sign convention for the potential 


V = —— + Vl m (51) 

r 

where the (1, m)-th spherical harmonic of the geopotential is given in plane polar coordi- 
nates by Kaula [9, eq. 3.61]: 

i 

V lm (r, it) = -pa* r _,_1 ^ F lmp Re(A; m exp i[(l - 2 p)u + m(£l - 0)]) (52) 

P= 0 


where u = u>+f and where again we have used the complex Ai m of Equation 6 to streamline 
the notation. In the plane polar coordinates the Lagrangian of an orbiting satellite is given 
by 


L = T-V= |r 2 + \r 2 ii 2 -V 

(53) 

from which we obtain the conjugate momenta 


dL . 
Pr = ~d> =T 

(54) 

dL 2 . 

(55) 

and the Hamiltonian 


H = T + V = \p 2 r + \r~ 2 p 2 u - pr _1 + V lm 

(56) 

Hamilton’s equations yield the equations of motion 


dH _ 3 2 _2 dVi m 

Pr= dr = r p u pr Qr 

(57) 

. dH 
r = = Pr 

dp T 

(58) 

, _ dH _ dv lm 

du du 

(59) 

. dH 2 

u = — = r *p u 

dp u 

(60) 

which are linearized to give the perturbation differential equations 


Ap r = (—2 r -3 - 3 r~Vu) Ar + (2r -3 p u ) A p u - 

(61) 

Ar = Ap r 

(62) 

.. dV, m 

Ap “ = du 

(63) 
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Aii = (—2r 3 p u ^ Ar + r 2 Ap u (64) 

Here the unperturbed or referenced system is determined by the Keplerian potential V = 
—pr~ l . The brackets “(•)” indicate time averaging, so that the bracketed expressions may 
be treated as constants. Equation 64 is superfluous. Since the perturbations are given in 
terms of the conjugate momenta instead of the velocities, the observation Equation 50 is 
restated as follows 

AR « sin — E(Ap r ) + r -1 cos —6(Ap u ) (65) 

We now compute the generalized disturbing force in (r, u)-coordinates. 

3V i 

= (k + 1 )p,a[r~ l ~ 2 £ F !mp Re(A im exp *[(/ - 2 p)u + m(Q - 0)]) (66) 

o r P = o 

£(* - 2 p)Fimp Irn (Aim exp i[(l - 2 p)u + m(Cl - 0)]) (67) 

p — 0 

In anticipation of performing a Fourier analysis, we use Hansen’s expansion to remove r 
and /, which are complicated functions of time, replacing them by a, which is constant, 
and M, which is linear in time. 

(-) exp i(l - 2 p)f = £ G ‘pi ex P - 2 P + <l) M (68) 

\a/ q 

(a) exp i(yl ~ 2p ^ = ^ G i+hP+\,g exp ~ 2p + ?) M ( 69 ) 

Note that for small e 

G,„(e) ~ ebl (70) 

The same holds for , of course. Substituting, we get: 

TTT = (* + £ F, mp Z G i ReS, m „ (71) 

or p=0 q '+ 1 .P+ 2 -? 

3V * 

= pala- 1 - 1 £(/ - 2p)Fim P £ G lpq Im S lmpq (72) 

OU p= 0 q 

where Si mpq is given by Equation 4. Note that u>, M, f2, and 6 are all linear in time for 
a Keplerian reference orbit. For polar orbits ft = 0. The “frozen orbit” condition lo = 0 
of Colombo [5, p. 117] would require some small expenditure of fuel. We now perform the 
Fourier analysis over a mission period, assuming an iVa-day and iV r -revolution stationary 
orbit. Let xp vary from 0 to 2tt over this period. 
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We now obtain some needed Fourier coefficient formulas. In addition to Equation 38, 
we have 

T [Re(w exp iwxf>) \ k] — + 6„, +fc u} (73) 

By requiring m ^ 0, / max < N r /2, and gcd(iV r ,iVd) = 1, the 6 w+ k terms vanish leaving 


and 


F[ImS lm M) | k] = -6 (l — 2 p+ q ) N t — mNd — k §1 mpq ( 0 ) 


T[ReSi mpq (^) | k] = -6(i — 2p + q ) N T — mNd — k *5 / mpq ( 0 ) 


These result in the Fourier coefficients of the generalized forces: 


T 


dV lm 
dr 1 

dV lm 


du 


» l 

= »a l e a 1 1 J2( l - 2 P) F ‘mrY, G ‘p<i^ S (i — 2p+ q ) N t — mNd — k ^Impq ( 0 ) 
p = 0 q M 


(74) 

(75) 


i i 

= (l + l) M a l e a-‘- 2 '£ F l m P E G l+1 + i (j ^(i— 2p+q)N T — mAT d — fc^/rap^O) (76) 

p=0 q 


We are now in a position to Fourier analyze the perturbation equations. 

ikipAp r (k) = (2 pr~ 3 - 3 r -4 />u) Ar(k) + (2r~ 3 p u ^ A p u (k) - ^y-(fc) 

ikxpAr(k) = Ap T (k) 
ikipAp u (k) = 

iktj>Au(k) = 2r -3 p u ^ A r(k) + r _2 A p u (k) 

Averaging the feedback coefficients, we can solve this system of algebraic equations for the 
spectra of the momentum perturbations. 


(77) 

(78) 

(79) 

(80) 
( 81 ) 


AMk) = t^ T+ (2pr-3 — 3r-rf) {< 2r ' 3p ") AP “ W “ TT^} (82) 

Ap *w - < 83 > 

This completes the consideration of the single satellite. The sum and difference identities 

^O _ O ( l ~ 2 P + <l) 8M O /OA 

Impq — £ COS 2 &lmpq (84) 

co ( l - 2 P + O / 0 c\ 

vSlmpq — sin ^ Slmpq (85) 
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are employed to obtain a final expression for the range rate spectrum: 


A R(k) = E fi(l-2p+q)N T -mN d -kFlmpSl mpq (0) 

Impq 

i _;_ 2 f . fiu (l — 2p + q)8M —ikxb 

X aa'a < sin — cos — : 

l 2 2 k 2 if>i + (2nr~ 3 - 3r- 4 pl) 


X pa e a < sin — cos 


(2 r 3 p u ) 


“ ( ' + 1)G '-h.»4.« 


, . 6u . ( l-2p + q)6Ml-2pG,„\ 

+ICOS _ sln — (86) 

This can be compared to the kinetic energy result Equation 43 by neglecting the terms 
proportional to sin and applying the substitutions: 


8u — 8M = iV r ^V> (88) 

= n/iV r (89) 

p u = a 2 n (90) 

A minor discrepancy remains, namely the factor of l — 2p appearing here in place of the 
factor / — 2p + q of the kinetic energy derivation. The consequences on observation matrix 
structure of using Equation 86 in place of Equation 43 are not severe, inasmuch as the 
non-zero elements are determined by the solution of the integer Equation 41 in either case. 
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3 Observation Equation 

3.1 Matrix Formulation 

The signal Equation 43 is linear in the unknown coefficients Ai m . This is easily seen by 
considering Equation 4. Note well that Ai m is merely a complex representation of the pair 
(C/ m , Sim ) according to Equations 6, 7, and 8; so that when we estimate Ai m , we will for all 
practical purposes be estimating C/ m and S; m . The signal equation thus becomes a linear 
observation: 

A R lm „(k) = + noise (91) 


We shall not at this time dwell on the nature of this noise term. For now we may consider it 
to be zero-mean gaussian with covariance a 2 J, where I is the identity matrix. Equation 91 
may be rewritten as 

A R(k) = £ dA f™ pq(k) A lm + noise (92) 

Impq * m 


where the summation indices are subject — at least in the case of m = 0 — to the 
constraint of Equation 41. Isolating the summations over l and m, we have: 

Afl(Jb) = ■£ dSm^Wfe) ^ + noiae (93) 

lm “Aim 

Since Equation 41 implies the mapping: 

k i— ► / — 2p + q,m (94) 

we can set . 

dA-R(fc) d XT pq ARlmpq(k) (95) 

dAlm dAlm K ’ 

where p and q are constrained by the fixed value of the expression l — 2p + q. In the 
special case of q = 0 considered below, we have p unique, making the summation over p 
and q a trivial one. The observation equation is now in matrix form and may be formally 
represented as 

y = Hx + v (96) 


Vk — ^ , Hh t i m Xi m -j- Vh 


where 


y k = A R(k) 

_ dAR(k ) 

— r\ A 

dA lm 

%lm = Al m 


— noise 
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3.2 Block Structure 


In this section we indicate how the observation matrix, i.e. dA R(k)/dAi m , decouples into 
independent blocks, and how these blocks factor into three matrix factors. Two of the 
factors, the row-dependent and column-dependent ones, are diagonal and the third (mid- 
dle) factor is real, in fact consisting of only Fi mp values. Finally, because of an additional 
symmetry, the number of rows of each block can be reduced by half. 

For simplicity, in this development we restrict the problem by requiring q = 0 and 
m ^ 0. Additional constraints are required by the method: 


gcd(IV r ,]V d ) = l 

(102) 

2 < l < N r /2 

(103) 

|*|<U*-(tfr + JVd) 

(104) 


Given these constraints, it has been shown in Thobe and Bose [14] that the linear dio- 
phantine equation 

k = (l — 2p)N T — mN<i (105) 

will have a unique solution for the unknowns l — 2 p and m, though such a solution may 
not be guaranteed to exist. A generalized Euclid’s algorithm is used in our software to 
solve Equation 105 for m and the set of l which satisfy it and the constraints. Since 
the rows of the observation matrix are indexed by k , this means that l — 2p and m are 
constant throughout any row. In fact, the matrix elements corresponding to a given order 
m and parity of degree |/[a = |/ — 2pj 2 form, after a suitable permutation of rows and 
columns, a block which is independent of all other orders and parities. Thus the solution 
of the observation equation may be carried out block by block, with massive computational 
savings. 

Given the above, an observation matrix element (see RHS of Equation 43) factors into 
three parts. 

= $row(fc)-Ff m p$col(0 (106) 

The first is dependent only upon the row index k. 


$row(fc) = cos • sin[Z - 2p]N T ^Y • - — exp ixpi mpO (0) (107) 

The second is simply the inclination function, which is real, and the the third is dependent 
only on the column multi-index (/, m). 



(108) 


Thus (each block of) the observation matrix is factored into the product of a diagonal 
matrix, a rectangular matrix, and another diagonal matrix. 
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Figure 1 shows the bounding non-zero elements of the observation matrix block m = 45, 
| / 1 2 == 0 . All even-even points on and inside the isosceles trapezoid are non- zero and all those 
outside are zero. For blocks such that |Z |2 = 1, the odd-odd elements are, of course, selected. 
The lowest and highest values of / are, respectively: 

min{/} = max{m, 2} + | max{m, 2} 4 - |/|2 | 2 (109) 

max{/} = /max - I /max + UU | 2 (HO) 


Since the oblique upper and lower boundaries of observation matrix block satisfy l—2p = 
±/, it is a simple matter to obtain the indices of the vertices. 

Having reduced the observation matrix essentially to a real isosceles trapezoid, we now 
explain the additional symmetry which allows the number of rows to be reduced by half. 
Consider what happens to the factors of Equation 106 as /— 2p — *• —(l—2p) or, equivalently, 
as p — > / — p. As can be seen from Figure 1, this constitutes a reflection of the matrix 
about the horizontal axis / — 2p = 0. We use the identity 

F, m[i _ p] (90°) = (— l)l'- 2 ^jr, mp (90°) (111) 


Defining k' = (— / + 2 p)N T — mNd and restricting / — 2p > 0, we have k representing rows 
in the top half of Figure 1 and k' representing their images in the bottom half. We then 
define the reduced measurement or observation vector to be: 


v(*) = 


1 

2 


$ row (*0 ^ v ' * row (fc') j 


Then solve the overdetermined linear system 


( 112 ) 


y(k) = 5Z Fi mp z(l, m) 

l,m 


(113) 


for z(l, m) and compute 

Ai m = z(l, m)/$coi(/, m ) (114) 

Note that the only significant computational effort is in solving Equation 113. This effort 
is an order of magnitude less than solving the original observation equation 


AR(k) = £ 

/,m 


dAR(k) 

dAi m 


Aim 


(115) 


which has twice the number of rows and is complex as well. 

The reduction of the observation equation just described can be rigorously proved 
using elementary properties of the pseudo-inverse. The reduced observation matrix for 
the block has the form of a right trapezoid. The computation of the reduced observation 
matrix elements Fi mp and the least squares solution of the reduced observation equation is 
discussed in Section 4. 
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Figure 1: Non-Zero Observation Matrix Elements for Block: m = 45, \l \2 = 0 
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4 Numerical Techniques 


4.1 Inclination Function Computation 

A stable, efficient, reliable method of computing the inclination functions Fi mp (i) is an 
essential and nontrivial element in the coefficient recovery process, as well as orbit predic- 
tion and simulation in a perturbed geopotential environment. We considered a number of 
existing techniques, all of which were unsatisfactory for our purposes for different reasons. 

Previously, we have not concerned ourselves with the question of normalization. We 
bring it up here because the use of the normalizing factor [9, formula (1.34)]: 


1 


( 7 + m)! 

(/ - m)!(2Z + 1)(2 - 6 m ) 


(116) 


sharply reduces numerical problems associated with excessive dynamic range in the val- 
ues of the harmonics and harmonic coefficients. While the other equations in this report 
are transparent with respect to whether normalized or unnormalized harmonics and har- 
monic equations are employed, actual expressions for the inclination functions Fi mp and 
the associated Legendre polynomials Pi m must be divided by Equation 116. The recovered 
harmonic coefficients will then be multiplied by the normalization factor. The expressions 
below are given in unnormalized form. The reader may apply the conversion himself. 

The first approach to computing Fi mp was by direct application of Kaula [9, formula 
(3.62)] 


(2/ - 2 ()! 


min{p,&} 

= g tltT- t)\(l - m - 2t)\2 2l ~ 2t 

min{p— £,/— m— 


sin 1 m 2t i 


(117) 


E 

5=0 


m 


cos t 


c=max{0,p— t— m-f s} 


( ! 


— m — 2t + s 
c 


)( p -- 3 c )(-'>'-* 


where k = [(/ — m)/ 2j (greatest integer or floor function). When specialized to t = 90°, 
the dummy index s can only have the value zero and the number of summations is reduced 
from three to two. 


min{p,fc} 


F lmr (W°) = ■£ 


(21 - 21 )! 


jz'q t\(l — <)!(/ — m — 2t)\2 2l ~ 2t 

x E 

c=max{0,p— t— m} 


(118) 


I 


m — 2t 

c 


)(,-? -c 


Even Equation 118 is computationally too slow for larger values of the indices /, m,p when 
large numbers of the coefficients are required, and we used it only for producing comparison 
values when testing our own method. 
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A second method which we attempted to use consists in part of evaluating P/ m at 
equally spaced intervals on an inclined circular orbit and employing an FFT to obtain the 
corresponding Fi mp . While interesting and accurate for small index values, say < 20; we 
found the implementation to be unstable for larger values, at least at the 90° inclination 
of interest. A further difficulty with the method is that the order in which the coefficients 
are produced: straight lexicographic ordering according to l,m,p (i.e. 200, 201, 202, 210, 
211, etc.). This would entail precomputing large numbers of coefficients, sorting them 
into blocks and rows according to the our observation matrix block structure, and storing 
them in a data base for recall during evaluation of the matrix elements. Our tests of this 
method were performed using a program kindly supplied by Clyde Goad of the Ohio State 
University Department of Geodetic Science. A description of this method is to be fo un d 
in Wagner [15]. 

To overcome these problems, we developed our own recursive method based on the an- 
alytical representation of Brumberg [3]. This new method is accurate, stable, and efficient 
to high index values. The inclination need not be restricted to 90°. Values are produced 
in a convenient order so that they can be computed on the fly, rather than precomputed, 
sorted, and stored. We have 


■ (sin ■ ,F, ( 

where 

s — 1 - (|Z - 2p - m| + |/ — 2p + m|)/2 

and 



i _ /_i \[^ 2i J+max{0,/-2p-m} (l)/+m( 2/ + 2p) max {(y_2p- TO }( 2p) max ^p i _/ + 2p+m} 

P ^ ^ ’ 2 , (l)p(l)|_ 2 p(l) | |_2p- m| 


(119) 

( 120 ) 


( 121 ) 


We derived Equation 121 ourselves from a formula of Izsak [8] according to Brumberg’s 
instructions when his own unnumbered formula preceding [3, eq. 22], gave incorrect values, 
probably due to a misprint. 

The notation (a) n = a(a + 1) ••• (a + n — 1) denotes the Pochhammer symbol or gen- 
eralized factorial, and 2 F 1 denotes the gaussian hypergeometric polynomial: 
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y (-a)n Q)n 3?” 

„0 Mn Ii! 


( 122 ) 


It turns out that direct application of Equation 122 leads to gross numerical difficulties, 
making it useless for all except calculations in a very small number of terms. If the eval- 
uation is carried out by simple evaluation and summing of the terms, then the method 
fails because the terms are very large, alternate in sign, and have a very small sum. If the 
evaluation is carried out by a simple recursion using Horner’s rule, then the method fails 
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because the resulting difference equation becomes unstable after a small number of itera- 
tions. According to Wagner [16], “fully recursive formulations in the literature . . . appear 
to be unstable at high degree.” 

Our contribution lies in inventing a recursive technique for evaluating 2 Fx, which is 
stable and efficiently produces accurate values for even large index values. Consider Equa- 
tion 119. The hypergeometric polynomial can be abbreviated as 


iF x 


— s, 21 — s + 1 
1 + \l — 2p — m\ 


sin 2 - = 2 Fi 


a,/3 

7 


x = F(a,/3) 


We apply two of the so-called gaussian contiguous relations 

(a - 0)F(a, 0) - aF(a + 1,0) + 0F(a, 0 + 1) = 0 


(123) 


(124) 


and 

(a - 0){1- x)F(a, 0) + (y- a)F(a; - 1, /?) - (7 - 0)F(a, 0- 1) = 0 (125) 

which may be found in Lebedev [11, p.242, eqs. 9.2.10-11]. In matrix form, these two 
second-order linear difference equations become: 


and 


" F(a,/? + 1) 


0—a a 
0 0 

F(a,/J) 

F(a„8) . 


1 0 

. F(a + 1,/?) 


(126) 


' F(a- 1,0 + 1) ' 


' fl^(l-z) 

1+1-7 I 

7-ar 

" F(a,/? + 1) 

F(a,0 + 1) 


1 

0 

F(a,/J) 


(127) 


Combining, we find that given F(a,0) and F(a + 1,0), we can by successively applying 
Equations 126 and 127, compute F(a — 1,0 + 1) and F(a, 0 + 1). Note that the transition 


(<*,$) ->(a- 1./3 + 1) 


(128) 


occurs when s — *• 5 + 1, which in turn occurs when l —t 1+2 and p — » p+ 1 simultaneously as 
m and l — 2p remain constant; in other words, as we take one step to the right along any row 
of the matrix block of Figure 1. Equations 126, 127, and 128 constitute a recursion, so that 
once starting values are computed using Equation 122, 2 Fi can be efficiently computed for 
an entire row. 

The best news is that this recursion turns out to be a stable one. This method has 
been used to evaluate the inclination function for index values up to one thousand. We 
used it to compute our observation matrix elements and also in our analytical trajectory 
simulation. 
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4.2 Observation Equation Solution 

Since the observation equation has been shown to have a block diagonal structure, it can 
be treated as a collection of uncoupled observation equations. Each block of the large 
observation matrix becomes the observation matrix for a smaller least squares estimation 
problem. Formally, let us write the (smaller) observation equation as 

y — Hx + v (129) 


For generality, regard all quantities as complex in the definitions: 

x = zero-mean n-vector of unknown coefficients 
y = zero-mean m-vector of measurements, rn > n 
v = zero-mean n-vector of measurement noises 
if = m x n observation matrix of rank m 

Since v is a unit -covariance complex white noise vector, it obeys the equations: 


Eu = 0 

(130) 

Era* = a 2 1 

(131) 

Eiw* = 0 

(132) 

That is to say, the real and imaginary parts of all components of v 

are mutually indepen- 


dent, and each has mean zero and variance <r 2 /2. 

The “best” solution of Equation 129 in the conditional mean sense, i.e. x = E{x | y}, 
is the one which conditionally minimizes the vector norm of the error: 

x = argmin{||y - Hx\\ 2 | y} (133) 


If U is a unitary (orthogonal in the real case) matrix, i.e. U 1 = 17*, then Equation 133 is 
equivalent to: 

(134) 


x = argmin |||y — Hx || 2 | yj 


where y = Uy and H = UH. Thus we have changed the problem but not the solution by 
multiplying both y and if by a unitary (or orthogonal) matrix U. 

To simplify bookkeeping, let us define the augmented matrix [if | y]. Then by judi- 
ciously selecting U, we can put [H \ y] into the form: 


U[H | y] = [UH | Uy] = [H\y) = 



(135) 
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where Hi is a full-rank upper triangular matrix and H^ is a zero matrix. 
From Equations 135 and 134 we have: 


x = argmin |||yi - #ia:|| 2 + ||y 2 || 2 | y l , y 2 } 

(136) 

Now Hi is invertible so that 


x = Hi l y i 

(137) 

since 


\\yi-Hixr + \M> = 0+\\y 2 \\ 2 

(138) 


which is clearly the minumum value. 

In practice, a finite sequence of unitary transformations is used: 


U = U (K) U (K ~ X) t/W (139) 

where 

tf( fc >[J?(*>|y( fc >] = [H^ly^] (140) 

For the Givens or rotation method, the elementary transformations are of the form: 


U (k) = 


a 


(k) 


a 


(k) 


a 


(k) 


a 


0 k ) 

33 


(141) 


where the ^-dependent indices i and j indicate the two rows of [H^ | affected by the 

transformation. Assuming that for some l > 1, the first l — 1 columns of are zero, and 
that for hffi and hffi are not both zero, then Equation 140 becomes: 


aW 

U u 

i 

8 


<*(*> 
L t 

a« 

33 J 



0 ••• 0 

0 ••• 0 


o Mr*’ ^ 

o >>f 


0 




/,(*+!) „(*+!) 1 

n in y% 

u( k + l ) 

L jn Uj 


( 142 ) 
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where 



Note that and h^ +1 ^ = 0. The reader may check for himself that this is true 

when for U ^ one selects the unitary matrix: 


a _ h^/rW hffi/rW 

a< $ a n - —hffi/rW h^/rW 


(144) 


The algorithm proceeds by selecting pairs of columns and rows to systematically zero 
out elements column by column from left to right, and row by row from the first subdi- 
agonal downward, until the pattern of Equation 135 is obtained. Then the solution x of 
Equation 137 is obtained by an efficient back substitution. So far we have described the 
standard method for solving least squares estimation problems by means of unitary Givens 
rotational transformations. For our problem, the augmented matrix has the special sparse- 
ness profile of a right trapezoid. For extra efficiency, one only need stop the downward 
scan when the bottom boundary of the profile is reached. It takes very little additional 
programming logic to accomplish this, which is one reason why the Givens method was 
selected. Furthermore, it is straightforward to include a priori estimates and covariances 
as well as correlated measurement noise. The main reason for selecting Givens rotations, 
however, is the excellent numerical stability of the method. 
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5 Simulation Results 


After having designed a spectral model for the satellite-to-satellite range rate, and after 
constructing software to recover the geopotential coefficients from the range rate data, it 
was necessary to exercise the software on simulation- produced data. That is the objective of 
this section. We will mainly concentrate on describing the tests and their results including 
mention of simplifying assumptions. Detailed descriptions of the math models and of the 
software algorithms are to be found in previous sections. 

The test article being subjected to validation testing is the recovery software itself. It 
is divided into parts according to the functions: 

• Obtaining the range rate spectrum from the data and classifying the spectrum into 
blocks according to order m and parity of degree |/| 2 ; 

• Generating the observation matrix for each block; 

• Constructing the observation vector and solving for the estimated geopotential har- 
monic coefficients. 

In order to exercise these functions and test the recovery software, we developed a one- 
sidereal day analytical simulation, restricting the maximum degree and order to eight. This 
allowed us to precisely control the test parameters so that discrepancies could be easily 
traced. Upon completion of this debugging and validation phase, we were able to proceed 
with confidence to tests involving foreign data, obtained from the much larger University 
of Texas simulation data set. Discrepancies at that stage could still be due to either our 
own modelling assumptions or to the Texas simulation, but not likely to errors in our 
software. We should point out that in accordance with good software engineering practice, 
detailed module testing of our software was performed prior to system level testing with 
any simulation data. 
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5.1 Analytical Validation 

A one day, 16-orbit trajectory of a satellite pair was simulated. The harmonic coefficients 
used are those of Rapp [12]. The perturbation included only tesseral terms, i.e. m ^ 0. As 
an additional check, the SST range rate was computed in two ways with identical results: 

• by a time domain computation followed by an FFT; 

• and by a direct frequency domain computation. 

The orbital parameters used in the simulation are given in Table 1. 


N r 

16 

number of revolutions 

N d 

1 

number of days 

a 

6 640 419m 

semi-major axis 

e 

0 

eccentricity 

i 

90° 

inclination 

u> 

7.292158 x 10“ 5 

perigee 

0 

90° 

ascending node 

Mo 

0 

initial mean anomaly 

R 

300km 

satellite separation 

a e 

6 378 155m 

earth radius 

n 

1.66745368 x 10" 3 

mean motion 


Table 1: Orbit Parameters for One-day Simulation 

The general time domain signal model equation is given by 

. • ... . ua l . _ 8M . ... „ 

ARlmpq(t) = -2 j^FlmpGlpg cos - 2 p+q)~z~] X 

w [f - 2p + q]M 
X [/ — 2p + q]M — mS ImS '”” (<) 


(145) 


Here we take M = N T tp and 9 = N d ip based on the assumption of a synchronized orbit and 
neglecting any secular motion. The frequency domain counterpart to Equation 145 with 
the noted assumptions is given by 


. • , IN . ual _ 8M . ... n . £M. 

^Rlmpqyk) — ® piq}~ I" 2 ^ ,rn P^*lpq COS — sill[(/ 2p <Jj “ ] X 


x [*-2p + g]M s (0) 

[/ — 2p -f q]M — mi) mp 9 


(146) 
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Figure 2: One-day SST Range Rate Perturbation Spectral Amplitude 

Summation in Equations 145 and 146 is performed over l, m, p, and q with the restrictions 
5 = 0 and k = (l — 2p + q)N T — mN a. 

Since this is the same model assumed by the recovery software, we expected the assumed 
coefficients to be recovered more or less exactly, which is precisely what occurred. As can 
be seen from Table 2, the difference between the model coefficients and the those recovered 
is seldom more than one or two in the fourth decimal place. The log amplitude (dB) of 
the magnitude of the spectrum of A R as a function of positive wave numbers k is given in 
Figure 2. 

Note that the structure of the spectrum is only clear when individual blocks are viewed. 
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Table 2: Recovered Coefficients for One-Day Analytically Simulated Mission 
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5.2 Texas Tape Recovery 

Our recovery software was exercised on SST range rate data produced at the University of 
Texas by Schutz et al. [13]. They simulated a 32-sidereal day mission of 525 near circular 
polar orbits. As required by our software, the orbit is synchronous — tuned to close on 
itself and repeat the same ground track each thirty two-day period. Table 3 gives the orbit 
characteristics. 


N r 

525 

number of revolutions 

N d 

32 

number of days 

a 

6 538 155m 

semi-major axis 

e 

0.00114 

eccentricity 

t 

90° 

inclination 

lj 

90° 

perigee 

n 

o 

O 

o 

ascending node 

M o 

0 

initial mean anomaly 

R 

300km 

satellite separation 

a e 

6 378 155m 

earth radius 

n 

1.1963697 x 10- 3 

mean motion 


Table 3: Orbit Parameters for University of Texas Simulation 

The data were supplied in ASCII in the form of a full 6250-bpi, 1800-ft magnetic 
tape. We converted this data to binary and extracted the SST range rate signal. The 
5-sec sampling interval was converted to 5.25903854 sec using linear interpolation so that 
precisely 2 19 = 524288 sampling intervals would equal one mission period. Due to the great 
length of the sample vector, we computed the complex spectrum of the range rate using 
the mass storage FFT of Fraser [7]. To assure ourselves of accurate numerical results, 
we carried out a simple analytical error analysis — converting the mass FFT software 
to double precision, comparing the results to those obtained with single precision, and 
applying the inverse FFT to obtain the original sampled signal. 

The spectrum was then permuted into blocks according to order m and parity of degree 
|/| 2 using a merge sort. Eack block was stored in a separate record of a direct-access file 
for easy retrieval. The recovery software computed each observation matrix block and 
solved the observation equation for that block thereby recovering the geopotential harmonic 
coefficients assumed by the University of Texas simulation. The equation solver had been 
tested on a 100 x 100 matrix of pseudo random numbers on the unit interval (0, 1]. When 
computing each block of the observation matrix, we also compute a predicted spectrum for 
that block using the assumed coefficients. This prediction should compare well with the 
spectrum extracted from the simulation data. We think it does, in particularly considering 
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that so far certain important considerations have been omitted from the analysis: 

• The whole range rate signal was used instead of the residual (perturbation) after 
removal of the reference orbit range rate. In effect an ideal circular reference orbit 
was used. 

• Secular perturbations of the orbit elements have not been included. This omission 
is justified by Kaula [9, eqs. 3.75, p. 40], in which the perturbation of Cl « 0 for an 
inclination of 90° and the perturbation of M + u> « 6.9° per day. 

• Terms in our signal model Equation 146 for which } / 0 have to date not been 
incorporated into the recovery software. 

Plots of the amplitude and phase of the SST range rate spectra are given in Fig- 
ures 26, 27, 5, and 6. The predicted spectra are those computed by us using the above 
described simplified kinetic energy model which is built into the recovery software. The 
observed spectra were computed by us from the University of Texas simulation data using 
the previously mentioned mass storage FFT and block permutation. 

It is interesting to note that the observed spectra degenerate into an approximately 
constant magnitude of — 130 dB and phase of ±60° beyond a wave number of approximately 
60 000 or 60% of bandwith of 100 260. We presume this reflects a limit to the fidelity of 
the U. of T. simulation. Our spectral model shows a definite and interesting structure in 
the outer 40% of the band, which is completely obscured in the simulation data. 

In the inner 60% of the band, the magnitudes agree fairly well, quantitatively as well 
as qualitatively. The main discrepancy in comparing both the magnitudes and phases is in 
the weakest points of the spectra. The phases agree best the closer one approaches k = 0. 
Out to about k = 20000 the discrepancy seldom exceeds 20°. 

It is likely that by including q = ±1 in the model, making improvements in the reference 
trajectory, considering secular motion, and iterating, that agreement would improve. The 
results so far are promising, indeed. 

Table 6 shows recovered values of the harmonic coefficients for m = 2 and / < 60 
and even. Except for l — 2, which would be heavily influenced by the uncorrected-for 
eccentricity in the reference orbit, the agreement between prediction and experiment is 
good. Also, no attempt was made to exclude the apparently erroneous data at extreme 
wave numbers. The computation was actually carried out to recover C/ m and Si m for 
values of l up to 180, but the higher l values were excluded for brevity. For values > 100 
the agreement is poor, since these values depend upon the extreme wave numbers. The 
software is general enough to be applied to any given block. The various refinements 
in modeling the signal spectrum would be expected to improve the agreement between 
predicted and estimated values. 
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6 Summary 

6.1 Conclusions 

We have described the analytical basis for our approach to the recovery of the harmonic 
coefficients of the geopotential field from satellite-to-satellite range rate data. We then pre- 
sented the results of numerical experiments with the objective of validating the approach, 
as well as demonstrating our ability to process the data on the scale required. 

The general basis of the method is that presented in Kaula [10]. However, the obser- 
vation equation has been freshly derived by us based upon much simpler expressions for 
the perturbation in the satellite kinetic energy, and the close relationship between varia- 
tions in kinetic energy and variations in satellite speed. This avoids much of the intricate 
dealings with perturbations of individual Keplerian elements. In addition, we borrowed 
from Colombo [5] the assumption of a synchronous orbit in which the ground track of each 
satellite repeats after N r revolutions and Nd days, these being a mutually prime pair of 
integers. This makes the spectrum discrete and allows one to employ fourier analysis to 
uncouple the observation equation into manageable-sized blocks. 

In the frequency domain, the observation matrix assumes a block diagonal structure. 
Without such a structure, the computational resources required to perform coefficient 
recovery would be truly massive. When q is restricted to zero as with circular polar orbits, 
there is one block for each combination of the harmonic order m and parity of degree 
|/| 2 . In addition, each block is further reduced to a real, right trapezoidal structure, which 
makes both its computation and solution even easier. 

On the numerical analysis front, we have developed a method, based on the mass- 
storage FFT of Fraser [7] and resampling of data, of obtaining the fourier coefficients of 
the sampled range rate over the entire mission period. A contribution which others in 
the orbit prediction field will find useful is the stable and efficient method of computing 
the inclination function Fi mp based on the analytical formulas of Brumberg [3] and the 
recursive generation of the hypergeometric polynomials using the contiguous relations of 
Gauss. We have employed the method of unitary Givens rotational transformations to 
solve the reduced and uncoupled frequency domain observation equations, while efficiently 
exploiting trapezoidal sparseness pattern of the blocks. 

We have demonstrated the low computational cost of the method. For a maximum 
degree of 180, each block can be set up and solved in a time less than or equal to approx- 
imately five minutes on a DEC VAX 11/750. While we have not yet reached the limits 
of accuracy of the recovered coefficients, we feel that the results achieved thus far on the 
University of Texas simulation data are highly encouraging. Experimental recovery of data 
generated by our own low degree and order analytical simulation shows excellent agree- 
ment between known and recovered coefficients, leading us to believe that considerable 
improvements in results of the Texas data experiments can be expected as well. 
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6.2 Recommendations 


It is recommended that the progress made this far be followed up and extended. A major 
priority is to improve the quality of the recovered coefficients. This can be done by: 

• improving spectrum modelling fidelity by: 

— employing the extended signal model of Section 2.3; 

— including the neglected q ^ 0 terms in the signal equation; 

— computing an eccentric reference orbit; 

— including secular perturbation effects. 

• employing iterative solution techniques including: 

— constructing intermediary reference orbits from modified Keplerian element per- 
turbations; 

— adding q ^ 0 terms to the observation matrix and employing operator pertur- 
bation methods; 

• extending preprocessing of recorded trajectory data using orbit determination tech- 
niques, filtering, and smoothing, to estimate mission period and other orbit charac- 
teristics. 

Our newly derived signal equation opens up the immediate possibility of investigating the 
following extensions of the method: 

• recovering the zonal (m = 0) coefficients; 

• using of common ground track trajectories rather than coplanar orbits; 

• using SST range rather than, or in conjunction with, the SST range rate as the 
physical measurement of choice. 

• performing parametric studies, e.g. varying the satellite-to-satellite separation angle; 

• tilting the orbit plane off of 90° inclination; 

• admitting non-zero eccentricities. 

We believe that additional simulation and validation will prove valuable. In particular we 
recommend: 

• perform our own Cowell method trajectory simulation; 

• use only small numbers of harmonics to limit simulation costs; 
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• use our efficient and stable computation of the inclination functions; 

• rule out possible errors in the University of Texas simulation data. 

In conclusion, we wish to stress the crucial role of effective data processing techniques in 
high resolution gravity mapping by satellite. 
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A Coefficients and Spectra 

This appendix contains 

• tables of known and recovered geopotential harmonic coefficients for degrees 2 through 
10 and orders 1 through 10; 

• magnitude plots of observed SST range rate spectra for degrees 2 through 10 and 
orders 1 through 10; 

• magnitude plots of predicted SST range rate spectra for degrees 2 through 10 and 
orders 1 through 10; 

• phase plots of observed and predicted SST range rate spectra for order 2 and even 
degree. 

The recovered data were obtained by processing the simulated range rate samples produced 
at the University of Texas by Schutz et al. [13]. 
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0.3783 xlO" 09 

-0.4727xl0-° 8 

0.4540 xlO" 09 

-0.4560x10"" 


B 

-0.1739 xlO" 08 

-0.2564 xlO" 08 

-0.1741 xlO" 08 

-0.3093x10-" 

El 

B 

0.8449 xlO' 08 

-0.7747 xlO" 08 

0.8799 xlO" 08 

-0.8503x10"" 

pa 

B 

-0.2872 xlO" 08 

-0.2495 xlO" 08 

-0.2939 xlO" 08 

-0.3018x10-" 

35 

1 


BsftlSfilsBlsBBi 


gnmiiEniBai 

37 

1 



0.5488 xlO" 09 

-0.1970x10-" 


B 

0.1054xl0-° 8 

0.4739 xlO" 08 

0.7538x10-" 

0.4863x10-" 


B 

-0.1244xl0-° 8 

-0.6224xl0-° 8 




B 

-0.5495 xlO" 09 

0.1426 xlO" 08 

0.1653xl0 _1 ° 

0.1368x10-" 

45 

1 

0.6479 xlO" 08 

-0.4578 xlO- 08 

0.6386 xlO" 08 

-0.4770x10-" 

47 

1 

-0.4427 xlO- 08 

-0.2029 xlO- 08 

-0.4816 xlO" 08 

-0.2427x10-" 

49 

1 

0.3104xl0-° 8 

0.1926xl0-° 9 

0.3324xl0 - ° 8 

0.5027 xlO -10 

51 

1 

0.1355xl0-° 8 

0.5262 xlO" 08 

0.1287xl0-° 8 

0.5867x10"" 

53 

1 

-0.1829xl0-° 9 

0.3123xl0-° 8 

-0.2536x10-" 

0.3811x10"" 

55 

1 

-0.3162 xlO" 08 

0.2393 xlO" 08 

-0.3048 xlO' 08 

0.2357x10-" 

57 

1 

0.3347 xlO" 08 

-0.2054xl0-° 8 

0.3466 xlO" 08 

-0.2877x10-" 

59 

1 



-0.3589 xlO- 08 

-0.9398 xlO' 09 

-0.4108 xlO" 08 

-0.1050x10-" 


Table 5: Recovered Coefficients for Block m = 1 , |/| 2 = 1 
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54 

2 

-0.3422 xlO" 09 

-0.9696 xlO" 10 

-0.1451 xlO" 08 

-0.2208 xlO- 09 

56 

2 

-0.4565 xlO" 08 

0.2627 xlO” 08 

-0.5034xl0-° 8 

0.3383 xlO" 08 

58 

2 

0.5005 xlO" 09 

0.2964xl0" 08 

0.8286 xlO" 09 

0.3406 xlO" 08 

60 

2 

0.3091 xlO" 08 

-0.1886 xlO" 08 

0.2667xl0 -08 

-0.2332 xlO" 08 


Table 6: Recovered Coefficients for Block m = 2 , 
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estimate 

truth 

l 

m 

C lm 

Sim 

Clm 

Sim 

3 

2 

■■imvmu 


■wnwriaitHa 


5 

2 





7 

2 



0.3664xl0 _O6 

0.1151 xlO" 06 

9 

2 

-0.3634xl0-° 8 

0.7449 xlO" 07 

0.3687 xlO" 07 

-0.3407 xlO" 07 

11 

2 

-0.1816xl0 _O7 

-0.1124xl0-° 7 

0.1644xl0-° 7 

-0.1000 xlO" 06 

13 

2 

0.2794 xlO- 08 

-0.5879 xlO" 08 

0.3141 xlO" 07 

-0.7884 xlO" 07 

15 

2 

-0.3252 xlO" 07 

0.4657xl0 _O7 

-0.9468 xlO- 08 

-0.1029 xlO" 07 

17 

2 







0.1326 xlO" 08 

0.4616 xlO" 07 

0.2100 xlO" 07 

0.1465 xlO” 08 



g>mur?f!niBa 




23 

2 

-0.1576 xlO' 07 

0.2733 xlO- 07 

-0.3115xl0 -09 

-0.9216 xlO- 08 

25 

2 




0.8387 xlO" 08 

27 

2 

-0.1168xl0-° 7 

0.3314xl0 -07 

0.1510xl0-° 8 

0.1559 xlO- 09 

29 

2 

-0.1302xl0-° 7 

■wmii 

-0.6303 xlO" 10 

-0.1084 xlO" 08 

31 

2 

-0.7414xl0-° 8 

0.3154xl0-° 7 

0.5693 xlO" 08 

0.2865 xlO- 08 

33 

2 

-0.1391 xlO- 07 

0.2816 xlO" 07 

-0.3023 xlO- 08 

0.1033 xlO' 08 

35 

2 

-0.2122 xlO' 07 

0.2803 xlO" 07 

-0.1125 xlO’ 07 

0.1704X10- 08 

37 

2 

-0.1022 xlO” 07 

0.1478 xlO" 07 

0.1144X10- 08 

-0.1125 xlO" 07 

39 

2 

-0.6168 xlO' 08 

0.3219 xlO” 07 

0.3519 xlO' 08 

0.7746 xlO" 08 

EM 

2 

-0.7113X10' 08 

0.2675 xlO" 07 

0.3365xl0 -08 

0.2446 xlO" 08 

ESI 








-0.1174xl0-° 7 

0.2107 xlO -07 

-0.2293 xlO' 08 


47 


-0.4169 xlO” 08 

0.2402 xlO“ 07 

0.5032 xlO" 08 


49 


-0.6067 xlO" 08 

0.2555 xlO" 07 

0.2439 xlO" 08 

0.3160 xlO’ 08 

El 

H 

-0.1340xl0-° 7 

0.1773xl0-° 7 

-0.3924 xlO" 08 

-0.3983 xlO" 08 


H 

-0.3291 xlO" 08 

0.2158xl0 _o7 

0.6056 xlO" 08 

-0.2101 xlO" 10 


B 


■inizkfiRinBa 

-0.2157X10- 08 

-0.3524 xlO' 08 


B 

-0.1125xl0- 07 

0.2149 xlO" 07 

-0.1878 xlO" 08 

-0.3267 xlO" 09 


El 

-0.5292 xlO' 08 

0.2250 xlO" 07 

0.4060 xlO" 08 

0.5905 xlO" 09 


Table 7: Recovered Coefficients for Block m = 2 , |/| 2 = 1 
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■ 


estimate 

truth 



Cl m 

Sim 

C lm 

Sim 

El 


0.1533 xlO- 06 

0.4716 xlO" 07 

0.9942 xlO’ 06 

-0.2022 xlO" 06 

Q 


Bt»KKKf?afigSl 


0.5620 xlO’ 07 

0.1750x10-°* 

8 

3 



-0.1624X10" 07 

-0.6824 xlO" 07 

10 

3 

-0.9960 xlO" 07 

-0.9940 xlO" 07 

-0.2184xl0-° 7 

-0.1500 xlO’ 06 

IEB 


0.7885 xlO" 09 

0.4635 xlO" 07 

0.5914xl0-° 7 

0.2832 xlO" 07 

Q 


0.4717x10-°* 

0.1526xl0-° 7 

0.4714xl0-° 7 

-0.4331x10-°* 

EH 


-0.5001 xlO" 07 

-0.1725x10-°* 

iHTTTHTm 

rnwH 

18 

B 

-0.2268 xlO" 07 

0.4357x10-°* 




B 

-0.2409 xlO" 07 

0.2070 xlO" 07 

-0.5467x10-°* 

0.1412X10" 07 



-0.6422 xlO" 08 

0.1650 xlO" 07 

0.8121x10-°* 

0.9310x10-°* 



-0.1897xl0-° 7 

0.5718x10-°* 

-0.6055x10-°* 

0.7817x10-°° 

EH 


-0.9953x10-°* 

0.8732x10-°* 

0.2038x10-°* 

0.3983x10-°* 



-0.9259x10-°* 

0.1307xl0-° 7 

WiitsBSfsaittoM 



Bl 

-0.9180x10-°* 

-0.1011 xlO" 07 

-0.8652x10-°° 

-0.1456 xlO" 07 


H 

-0.1255 xlO- 07 

6.5310x10-°* 

-0.5179x10-°* 

0.2764x10-°* 

EH 


0.4211x10-°* 

0.9458x10-°* 

0.1227xl0-° 7 

0.6631x10-°* 

36 

3 

-0.4421x10-°* 

-0.6014x10-°* 

0.1311x10-°* 

-6.9415x10-°* 



-0.1700x10-°* 

0.3068x10"°* 

0.4206x10"°* 

0.7150x10-°° 



-0.5889x10-°* 

0.3516x10-°* 



42 

3 

-0.3014x10-°* 

0.1044xl0-° 7 


1 

44 

3 

-0.7043xl0 _1 ° 

-0.6512x10"°* 

0.4865x10“°* 

-0.9264x10-°* 

Em 


-0.5420x10-°* 

0.1444x10-°* 

-0.9530x10-°° 

0.1723x10-°° 



-0.4603x10-°* 

0.7096x10-°° 

-0.7420x10"°° 

-0.1397x10-°* 

50 

3 

-0.4063x10-°* 

0.2287x10-°* 

-0.9008 xlO" 10 

-0.6754xl0 -1 ° 

52 

3 

-0.2204x10-°* 

-0.1738x10-°* 

0.2020x10"°* 

-0.2870x10-°* 

54 

3 

-0.3967x10-°° 

0.3214x10-°* 

0.4001x10-°* 

0.1551x10-°* 

56 

3 

-0.5718x10-°* 

0.3769x10-°* 

-0.3507x10-°* 

0.2631x10-°* 

58 

3 

-0.6050x10-°* 

-0.2390x10-°* 

-0.2602x10-°* 

-0.5124x10-°* 

60 

3 

0.4281x10-°° 

0.1946x10"°* 

0.3904x10-°* 

0.1673x10-°* 


Table 8: Recovered Coefficients for Block m = 3 , |Z |2 = 0 
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estimate 

truth 

l 

m 

c lm 

Sim 

Clm 

Sim 



0.6979 xlO" 06 

0.1167xl0-° 5 

0.7186 xlO" 06 

0.1417xl0-° 5 



-0.2431 xlO- 06 

-0.4817xl0-° 7 

-0.4561 xlO" 06 

-0.2175 xlO" 06 



0.2714xl0 _O6 

-0.1019 xlO' 06 

0.2413X10- 06 

-0.2066 xlO- 06 

9 

3 

-0.1122 xlO" 06 

0.1455 xlO" 08 

-0.1696 xlO" 06 

-0.4890 xlO- 07 

11 

3 

-0.2274 xlO" 07 

-0.8263 xlO" 07 

-0.5211 xlO- 07 

-0.1181 xlO- 06 

m 


-0.5361 xlO" 08 

0.8903 xlO' 07 

-0.2617xl0-° 7 

0.7483 xlO" 07 

[15 

3 

0.5755 xlO" 07 

0.3770 xlO- 07 

0.4578 xlO" 07 

0.2490 xlO" 07 

Bfl 

D 

0.1247 xlO" 07 

0.1914xl0-° 7 

0.3229 xlO- 08 

0.9191 xlO- 08 

EES 

■a 

-0.3787 xlO -08 

-0.2998 xlO" 08 

-0.1299 xlO’ 07 

-0.1209 xlO" 07 

El 


0.2096 xlO" 07 

0.2011 xlO" 07 

0.1550 xlO" 07 

0.1320 xlO" 07 

EM 


-0.7860 xlO" 08 

-0.7228 xlO" 08 

-0.1499 xlO” 07 

-0.1462 xlO” 07 



-0.8604xl0 - ° 8 

0.5417xl0-° 8 

-0.1387xl0-° 7 

0.3318 xlO" 09 

gn 


0.6637xl0-° 8 

0.6051 xlO- 08 

0.2441 xlO- 08 

0.1593 xlO" 08 



0.4905 xlO" 08 

-0.1179 xlO" 08 

0.1166xl0-° 8 

-0.6380 xlO" 08 

31 

3 

BsEsUsBIsBIH 


-0.3106 xlO- 08 

-0.8561 xlO" 08 

Mm 


-0.1281 XlO- 08 

0.5220 xlO" 08 

-0.3819 xlO' 08 

0.2197 xlO" 08 


El 

0.9278 xlO" 08 

0.3788 xlO- 08 

0.6000 xlO" 08 

0.1159xl0-° 9 

m 

B 

0.2585 xlO” 08 

0.2392 xlO’ 08 

0.8462xl0- 10 

0.5645 xlO" 09 

Kil 

D 

0.5326 xlO" 09 

0.8049 xlO' 08 

-0.2278 xlO" 08 

0.5925 xlO- 08 

41 

B 


gilifeHifilftBSa 





0.1228xl0-° 8 

0.1442 xlO- 08 

-0.8106 xlO" 09 

-0.1860 xlO- 08 

45 

3 



-0.2299 xlO- 08 

-0.4414 xlO” 08 

47 

3 



0.7248 xlO" 09 

0.2860 xlO- 08 

49 

3 



-0.2684 xlO" 08 

0.4022 xlO’ 08 

EM 




-0.4379 xlO- 08 

-0.6072 xlO" 08 

El 


-0.7821 xlO" 09 

0.3073 xlO" 08 





0.6751 xlO" 08 

0.2234 xlO" 08 

0.5257 xlO" 08 

0.6296 xlO" 09 



-0.2258 xlO" 08 

0.5332 xlO" 08 

-0.3694 xlO" 08 

6.3792 xlO" 08 

59 

3 



0.5700 xlO- 09 

-0.5525 xlO- 08 


Table 9: Recovered Coefficients for Block m = 3 , |/| 2 = 1 


42 
















estimate 

truth 

l 

m 

c, m 

Sim 

Clm 

Sim 

4 

4 

rimmsa 

-0.2678 xlO" 06 



6 

4 

-0.8381 xlO" 08 

-0.4253 xlO' 06 

-0.9175 xlO- 07 

-0.4687 xlO" 06 

8 

4 

-0.1603 xlO- 06 

0.5026 xlO” 07 

-0.2438 xlO" 06 

0.7621 xlO" 07 

m 

4 

-0.5495 xlO~ 07 

-0.8789 xlO" 07 

eniiHKEinBia 








14 

4 

0.1721 xlO” 07 

-0.2953 xlO" 07 

0.3440 xlO" 08 

-0.2587xlO- 07 

m 


■tEKwranga 


0.3724xl0-° 7 

0.6272 xlO- 07 

18 

4 

| 


wminvim 


20 

4 

0.7858 xlO' 08 

-0.2046 xlO" 07 

0.2339 xlO” 08 

-0.1833 xlO” 07 





0.6285 xlO" 08 

0.3154xl0-° 7 

24 

4 

0.8392 xlO" 08 

0.2735 xlO" 08 

0.6034 xlO- 08 

0.4618 xlO" 08 





0.9073 xlO" 08 

-0.5606 xlO" 08 

28 

4 

0.4827 xlO" 08 

0.3404 xlO” 08 

0.2244xl0-° 8 

0.4970 xlO” 08 

30 

4 

-0.1563 xlO” 08 

-0.3982 xlO" 09 

-0.1821 xlO" 08 

0.8802 xlO" 09 

32 

4 

0.2723 xlO” 08 

-0.1024xl0-° 7 

0.3081 xlO" 09 

-0.9114xl0-° 8 

34 

4 

-0.8392 xlO- 09 

-0.1329 xlO" 08 

-0.2465 xlO- 08 

-0.2254xl0-° 9 

36 

4 

0.2712 xlO” 08 

-0.1896 xlO" 08 

0.1652xl0-° 8 

-0.1150 xlO” 08 

38 

4 

0.4448 xlO- 08 

-0.1371 xlO- 08 

0.2565 xlO" 08 

-0.5651 xlO" 09 

40 

4 

0.2511 xlO' 08 

-0.8719xl0-° 8 

0.1973 xlO- 08 

-0.7505 xlO” 08 

42 

4 

0.2387 xlO" 08 

0.8158 xlO" 09 

0.1846xl0-° 8 

0.1485 xlO" 08 

44 

4 

0.4612 xlO" 08 

-0.8832 xlO" 09 

0.3702 xlO” 08 

-0.2616 xlO" 09 

46 

4 

0.2954 xlO" 08 

-0.6171 xlO" 08 

0.1743 xlO" 08 

-0.5470 xlO' 08 

48 

4 

-0.7957 xlO" 09 

-0.3980 xlO" 09 

-0.1717xl0-° 8 

0.2607 xlO" 09 

50 

4 

-0.1026 xlO" 07 

0.2690 xlO' 08 

-0.1191 xlO" 07 

0.2552 xlO" 08 

52 

4 

0.3199 xlO" 08 

0.1057 xlO" 08 

0.2909 xlO" 08 

0.1443 xlO" 08 

54 

4 

0.5895 xlO" 09 

-0.2972 xlO" 08 

-0.6773 xlO" 09 

-0.2658 xlO" 08 

56 

4 

0.8456 xlO" 09 

0.3098 xlO" 08 

-0.8299 xlO" 10 

0.3959 xlO' 08 

58 

4 

-0.1509 xlO" 08 

-0.2837 xlO" 08 

-0.1607xl0 -08 

-0.1826 xlO" 08 

60 

4 

0.7255 xlO" 08 

0.1546 xlO" 08 

0.7021 xlO" 08 

0.2241 xlO" 08 


Table 10: Recovered Coefficients for Block m = 4 , |/| 2 = 0 
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estimate 

truth 

l 

m 

Cl m 

Sim 

Clm 

Sim 

5 

4 



-0.2899 xlO” 06 

0.4979 xlO" 07 

7 

4 

-0.2488 xlO- 06 

0.1012 XlO- 07 

-0.2721 xlO- 06 

-0.1281 xlO” 06 

9 

4 

fmmai 

■nmrftEinBa 

-0.1425 xlO” 07 

0.3122xl0-° 7 

11 

4 

-0.7311 xlO" 07 

-0.1180xl0-° 7 

-0.4948 xlO" 07 

-0.9853 xlO" 07 



-0.3162 xlO” 07 

0.4863 xlO" 07 

-0.3766 xlO" 08 

-0.1850 xlO” 07 

15 

4 

-0.5069 xlO' 07 

0.5675 xlO" 07 



17 

4 

-0.3419 xlO" 07 

0.6640 xlO" 07 



19 

4 

-0.9324X10" 08 

0.3870 xlO” 07 



21 

4 

-0.2108 xlO" 07 

0.4319 xlO" 07 







-0.1119xl0-° 7 

0.6884xl0-° 8 

25 

4 

-0.1704xl0 -07 

0.3248 xlO" 07 

0.7004xl0 -09 

-0.1541 xlO’ 08 

27 

4 







-0.3953 xlO- 07 

0.3034 xlO" 07 

-0.2486 xlO- 07 

0.3877 xlO -09 

m 



■n 

0.8584xl0-° 8 

-0.5909 xlO” 08 

33 

4 

-0.1676 xlO- 07 

0.2725 xlO" 07 

-0.1881 xlO- 08 

0.8881 xlO” 09 




ITO'ffrltf 

-0.1558 xlO" 08 

0.1530 xlO" 08 




wvwmm 

0.5658 xlO" 08 

-0.1382xl0-° 8 



-0.2107 xlO” 07 

0.1652 xlO- 07 

-0.7801 xlO" 08 

-0.8662 xlO" 08 

41 

4 

-0.1737 xlO" 07 

0.2564xl0-° 7 

-0.4081 xlO- 08 

0.2775 xlO’ 08 

43 

4 

-O.llllxlO -07 

0.2383 xlO" 07 

0.2402 xlO" 08 

0.3607 xlO- 09 

45 

4 

-0.9467 xlO" 08 

0.2036 xlO" 07 

0.3647 xlO" 08 

-0.2942 xlO" 08 

47 

4 

-0.1293 xlO" 07 

0.2511 xlO" 07 

-0.1741 xlO" 09 

0.3373 xlO" 08 

49 

4 

-0.1374X10- 07 

0.3068 xlO" 07 

-0.1172xl0-° 8 

0.9215 xlO" 08 

51 

4 

-0.1406 xlO" 07 

0.2433 xlO' 07 

0.1636xl0-° 9 

0.2250 xlO" 08 

53 

4 

-0.9025 xlO" 08 

0.2112xl0-° 7 

0.3625 xlO" 08 

-0.9974xl0-° 9 

55 

4 

-0.1262 xlO" 07 

0.2008 xlO' 07 

0.1129xl0-° 9 

-0.1052 xlO" 08 

57 

4 

-0.1736 xlO" 07 

0.1678 xlO" 07 

-0.5386 xlO" 08 

-0.5801 xlO" 08 

59 

4 

-0.1013 xlO" 07 

0.2140 xlO" 07 

0.3325 xlO" 08 

0.4438 xlO' 09 


Table 11: Recovered Coefficients for Block m = 4 , |/|2 = 1 
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/ 

rrt 

6 

5 

8 

5 

10 

5 

12 

5 

14 

5 

16 

5 

18 

5 

20 

5 

22 

5 

24 

5 

26 

5 

28 

5 

30 

5 

32 

5 

34 

5 

36 

5 

38 

5 

40 

5 

42 

5 

44 

5 

46 

5 

48 

5 

50 

5 


estimate 


truth 


54 5 


0.3213 xlO- 07 


0.9594xl0-° 7 

0.6149xl0 -07 


0.1843 xlO" 07 
0.2549 xlO" 07 
0.1688xl0~° 7 


0.1204xl0 _O7 
0.8765 xlO" 08 


0.9437xl0-° 7 


0.9601 xlO' 07 
0.5106 xlO- 07 


0.5881 xlO" 07 
0.4825 xlO" 07 
0.1947xl0-° 7 


0.2550 xlO* 07 
0.1440 xlO" 07 


0.8637xlO -08 0.6605 xlO' 08 


0.4313 xlO" 08 0.1409 xlO" 07 


0.9098 xlO" 08 0.6334 xlO" 08 


8 


lEESSEliliiSi 

■ 


0.5624xl0-° 8 0.2034xl0"° 8 


8 


-0.2638 xlO” 06 -0.5362 xlO” 06 


-0.1945 x: 
-0.5876 xlO" 07 
0.4510 xlO" 07 
0.2296 xlO" 07 
-0.1170x: 
0.3841 x: 


.8329 x: 


-0.4115xl0-° 7 
0.8910xl0-° 8 
-0.1503 xlO" 07 


-0.3324xl0 _O8 0.7321 xlO- 09 


-0.4508 xlO' 08 -0.8155 xlO- 08 
0.7388 xlO" 08 0.3360 xlO- 08 






0.1503xl0-° 8 -0.2284xl0 _O8 



45 
































■ 


estimate 

truth 

B 


Cl m 

Sim 

Clm 

Sim 


□ 

Kmmmivssm 


0.1733 xlO" 06 



D 



0.1253 xlO" 07 


9 

5 

0.4263 xKT 07 

-0.3675 xlO" 07 

-0.2726 xlO- 07 

-0.3387 xlO" 07 

11 

5 

0.7550 xlO" 07 

0.1392 xlO" 07 

0.4477 xlO" 07 

0.3494 xlO" 07 

13 

5 

0.8485 xlO" 07 

0.5233 xlO- 07 

0.6588 xlO" 07 

0.7450 xlO' 07 

15 

5 

0.2214xl0 -07 

-0.8543 xlO" 08 

0.4212 xlO" 08 

0.6387 xlO" 08 

17 

5 


Emma 

-0.1595 xlO~ 07 

0.8627 xlO" 08 

19 

5 

0.1755 xlO" 07 

-0.1232 xlO" 08 

0.7192 xlO" 08 

0.8381 xlO" 08 

21 

5 

0.1159 xlO" 07 

-0.1488 xlO" 07 

0.3017 xlO" 08 

-0.6414 xlO- 08 

23 

5 

0.8889 xlO" 08 

0.1707 xlO" 08 

0.1813 xlO" 08 

0.9295 xlO- 08 

25 

5 

bjsbBKsEI 


-0.7896 xlO -08 

-0.2705 xlO" 08 

27 

5 

0.2195 xlO" 07 

0.3084xl0" 08 

0.1721 xlO' 07 

0.8708 xlO" 08 

29 

5 




nvgnQBB|| 

31 

5 

0.2946 xlO* 08 

-0.3015 xlO" 08 

-0.1341 xlO" 08 

0.2176 xlO" 08 

33 

5 

0.3800 xlO“ 08 

0.1326 xlO" 08 

-0.3341 xlO" 09 

0.5819 xlO' 08 

35 

5 

-0.7749 xlO' 09 

-0.6035 xlO' 08 

-0.4412 xlO" 08 

-0.2658 xlO" 08 

37 

5 

-0.2644 xlO" 08 

0.3774xl0-° 8 

-0.6016 xlO" 08 

0.7904 xlO" 08 

39 

5 



0.6153 xlO" 08 

0.7278 xlO" 08 

41 

5 



HrrmriB 

KgEBfQnsa 

43 

5 

-0.9821 xlO' 08 

-0.1264 xlO" 08 

-0.1252 xlO“ 07 

0.2722 xlO* 08 

45 

5 

0.5496 XlO" 08 

-0.2271 xlO" 08 

0.2959 xlO' 08 

-0.1212 XlO" 09 

47 

5 

0.2656 xlO' 09 

-0.5855 xlO" 08 

tHTrTfflll 


49 

5 





51 

5 





53 

5 

0.1896 xlO" 08 

-0.5368 xlO" 08 

-0.9796 xlO" 10 

-0.3444 xlO" 08 

55 

5 

0.9494 xlO" 08 

0.6484 xlO" 08 

0.8351 xlO" 08 

0.9581 xlO" 08 

57 

5 

-0.1994xl0~ O8 

0.3272 xlO" 08 

-0.4272 xlO' 08 

0.5538 xlO" 08 

59 

5 

0.3510 xlO- 08 

-0.1428 xlO- 08 

0.1513 xlO" 08 

0.9107 xlO" 09 


Table 13: Recovered Coefficients for Block m — 5 , |Z |2 = 1 
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estimate 

truth 

l 

m 

Clm 

Sim 

Clm 

Sim 

6 

6 

0.1750 xlO" 06 

-0.7622 xlO' 07 

0.6872 xlO' 08 

-0.2371 xlO" 06 

8 

6 

0.2833 xlO" 07 

0.2032 xlO" 06 

-0.6642 xlO" 07 

0.3110xl0-° 6 

10 

6 

0.1275xl0 -07 

-0.1049 xlO" 06 

-0.4372 xlO' 07 

-0.8351 xlO’ 07 

12 

6 

0.3468 xlO" 07 

0.1090 xlO" 07 

0.5968 xlO" 08 

0.3544 xlO" 07 

14 

6 

0.5246 xlO" 08 

-0.1000 xlO" 07 

-0.1341 xlO" 07 

0.4073 xlO" 08 

16 

6 

0.1428 xlO" 07 

-0.5287xl0-° 7 

- 0.1137xl0-° 8 

-0.4556 xlO" 07 

18 

6 

0.2102 xlO" 07 

-0.2510 xlO' 07 

0.1421 xlO" 07 

-0.1700 xlO" 07 

20 

6 

0.1826 xlO" 07 

-0.6249 xlO" 08 

0.1317xl0-° 7 

-0.9172 xlO" 09 

22 

6 

0.1649 xlO" 07 

-0.7016 xlO" 08 

0.1266 xlO' 07 

-0.1379 xlO" 08 

24 

6 

0.4380 xlO" 08 

-0.3599 xlO" 08 

O.HOOxlO" 08 

0.4751 xlO" 11 

26 

6 

O.H72xlO -07 

-0.7606 xlO" 08 

0.8964xl0-° 8 

-0.3752 xlO" 08 

28 

6 

0.6150 xlO" 08 

0.1314xl0-° 8 

0.3591 xlO' 08 

0.5057 xlO" 08 


6 

0.5871 xlO" 08 

0.2581 xlO" 08 

0.3355 xlO" 08 

0.6076 xlO- 08 

32 

6 

0.1741 xlO" 08 

-0.7688 xlO" 08 

-0.1762 xlO" 09 

-0.5283 xlO” 08 

34 

6 

0.6993 xlO" 08 

0.4661 xlO" 08 

0.5561 xlO" 08 

0.8176 xlO” 08 

36 

6 

0.1488 xlO" 07 

-0.3454 xlO" 08 

0.1383 xlO’ 07 

-0.1517xl0-° 8 

38 

6 

-0.1142 xlO" 07 

0.1617xl0-° 8 

-0.1371 xlO" 07 

0.3657 xlO" 08 



-0.5647 xlO" 09 

0.1129xl0-° 8 

-0.1484xl0-° 8 

0.3410xl0-° 8 

42 

6 

0.3408 xlO" 08 

-0.4757 xlO” 08 

0.2664 xlO' 08 

-0.3197xl0-° 8 

44 

6 

BnwteRgiiiBBa 


-0.7911 xlO" 08 

0.2485 xlO” 08 



wmmm 


-0.4882 xlO" 08 

-0.1506 xlO" 08 





0.4229 xlO" 08 

0.5803 xlO" 08 

50 

6 

■iBEQGHlSil 


0.4259 xlO" 09 

0.1113xl0-° 8 



-0.4011 xlO" 08 

-0.3468 xlO" 08 

-0.5492 xlO" 08 

-0.1997xl0-° 8 

54 

6 

0.4193 xlO" 08 


0.3399 xlO" 08 

-0.2739 xlO' 08 

56 

6 

-0.4150 xlO" 08 

-0.5174xl0-° 9 



58 

6 

0.9002 xlO" 10 

-0.1372xl0-° 8 

-0.1102xl0-° 8 

-0.2030 xlO" 10 

60 

6 

-0.4797 xlO" 08 

-0.1699 xlO" 08 

-0.5000 xlO" 08 

-0.3095 xlO" 09 


Table 14: Recovered Coefficients for Block m = 6 , |/| 2 = 0 
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estimate 

truth 



Clm 

Sim 

Clm 

Sim 

7 

6 

-0.1294xl0 _O6 

0.2017 xlO' 06 

-0.3607X10' 06 

0.1518xl0'° 6 

9 

6 

0.9954xl0'° 7 

0.2821 xlO' 06 

0.5770 xlO' 07 

0.2106 xlO' 06 





-0.9824 xlO' 08 

0.3571 xlO' 07 

13 

6 

-0.1714X10- 07 

0.1070xl0'° 6 

-0.2349 xlO' 07 

0.2598 xlO' 08 

15 

6 

0.1979 xlO" 07 

0.6086 xlO' 07 

0.2448 xlO' 07 

-0.4019 xlO' 07 

17 

6 

-0.2125 xlO' 07 

0.5479 xlO' 07 

-0.1682xl0'° 7 

-0.3550 xlO' 07 

19 

6 

-0.1073 xlO" 07 

0.1011 xlO' 06 

-0.3979 xlO' 08 

0.2140 xlO' 07 

21 

6 

-O.lllOxlO" 07 

0.7648 xlO' 07 

-0.3420 xlO' 08 

0.7638 xlO' 09 

23 

6 

-0.1717X10" 07 

0.9139 xlO' 07 

-0.9506 xlO' 08 

0.2215 xlO' 07 

25 

6 

-0.2945 xlO" 08 

0.6887 xlO' 07 

0.6582 xlO' 08 

0.2377 xlO' 08 

27 

6 

-0.3388 xlO" 08 

0.6699 xlO' 07 

0.5548 xlO' 08 

0.4585 xlO' 08 

29 

6 

-0.2295 xlO" 08 

0.6364 xlO' 07 

0.7066 xlO' 08 

0.3580 xlO' 08 

31 

6 

-0.1018xl0"° 7 

0.6043 xlO' 07 

-0.1459 xlO' 08 

0.2522 xlO' 08 

33 

6 

-0.9110xl0'° 8 

0.5340 xlO' 07 

-0.1414xl0' n 

-0.1530 xlO' 08 

35 

6 

-0.5891 xlO" 08 

0.6033 xlO' 07 

0.3154xl0'° 8 

0.6401 xlO' 08 

37 

6 

-0.9655 xlO' 08 

0.5766 xlO' 07 

-0.2523 xlO' 08 

0.6194xl0'° 8 

39 

6 

-0.1432 xlO' 07 

0.5178xl0~° 7 

-0.4754 xlO' 08 

0.8801 xlO' 09 

41 

6 

-0.8327xl0' 08 

0.5094 xlO' 07 

0.5032 xlO' 09 

0.1017xl0'° 8 

43 

6 

0.5052 xlO' 09 

0.4719xl0'° 7 

0.9180 xlO' 08 

-0.1126xl0' 08 

45 

6 

-0.1122 xlO' 07 

0.4578 xlO' 07 

-0.2392 xlO' 08 

-0.2470 xlO' 08 

47 

6 

-0.6049 xlO' 08 

0.4613 xlO' 07 

0.3306 xlO' 08 

-0.7040 xlO' 09 

49 

6 

-0.7116xl0' 08 

0.4899 xlO' 07 

0.8487xl0'° 9 

0.2003 xlO' 08 

51 

6 

-0.6354xl0' 08 

0.4504 xlO' 07 

0.1865 xlO' 08 

-0.1430xl0'° 8 

53 

6 

-0.5773 xlO' 08 

0.4538 xlO' 07 

0.3644xl0'° 8 

-0.2797 xlO' 09 

55 

6 

-0.6518 xlO' 08 

0.4507 xlO' 07 

0.1263 xlO' 08 

-0.5922 xlO' 10 

57 

6 

-0.7728 xlO' 08 

0.4444 xlO' 07 

0.1513 xlO' 08 

-0.1098 xlO' 08 

59 

6 

-0.1423 xlO' 07 

0.4395 xlO' 07 

-0.6645 xlO' 08 

-0.1370 xlO' 08 


Table 15: Recovered Coefficients for Block m = 6 , |/| 2 = 1 
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estimate 

truth 

l 

m 

c lm 

Sim 

Clm 

Sim 

8 

7 


EjESQESH 



10 

Kfl 

-0.5507 xlO- 07 

-0.8564xl0- 07 




n 



-0.2024xl0- 07 

0.2704 xlO' 07 

14 

7 

■nan 



■WlffillH 

16 

Kfl 



-0.9166 xlO" 08 

-0.1529 xlO" 07 


m 



■iiwawgntHa 



7 





22 

7 





24 

7 



m — 


26 

7 





28 

7 

-0.9597 xlO" 08 

-0.7304 xlO” 08 



30 

7 

-0.7250x10"" 

-0.8703 xlO- 08 

0.7415 xlO* 08 

0.6068x10-" 

32 

7 


rTrmrlaa 


■rmitriii 

34 

7 



-0.3106 xlO" 09 

-0.1337x10-" 

36 

7 

-0.6353 xlO" 08 

-0.1048 xlO" 08 

-0.4350x10-" 

0.5864x10"" 

38 

7 



-0.1142 xlO- 08 

-0.8848x10-" 

40 

7 

-0.6942 xlO" 08 

-0.7329 xlO" 08 

-0.2567 xlO* 08 

-0.1286x10"" 

42 

7 

-0.1611 xlO" 08 

-0.1123 xlO" 07 

0.2737 xlO" 08 

-0.6817x10"" 

44 

7 

0.2337 xlO" 08 

0.3946 xlO' 08 

0.7088 xlO" 08 

0.9518x10-" 

46 

7 

-0.2386 xl0~" 

-0.1580 xlO" 07 

0.3778 xlO" 08 

-0.1278 xlO" 07 

48 

7 

-0.6055 xlO" 08 

-0.4999 xlO" 08 

-0.2541 xlO' 08 

0.7060 xlO" 10 

50 

7 

-0.6774 xlO" 09 

-0.1839 xlO" 09 

0.3055xl0-° 8 

0.4305x10"" 

52 

7 

-0.6316 xlO" 08 

-0.4045 xlO" 08 

-0.3350x10-" 

-0.1222x10-" 

54 

7 

-0.1410 xlO" 08 

-0.6386 xlO" 08 

0.1933x10-" 

-0.1875x10-" 

56 

7 

-0.1567xl0-° 8 

-0.3518 xlO" 08 

0.1196x10-" 

-0.3526x10"" 

58 

7 

-0.8121 xlO’ 08 

0.8855 xlO" 09 

-0.5674x10-" 

0.4590x10"" 

60 

7 

-0.2759 xlO" 08 

-0.4326 xlO- 08 

0.4289x10-" 

-0.1701x10-" 


Table 16: Recovered Coefficients for Block m = 7 , [Z| 2 = 0 
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Table 17: Recovered Coefficients for Block m = 7 , |/| 2 = 1 


50 























estimate 

truth 



C tm 

Sim 

C/m 

Sim 

8 

8 

0.2016 xlO- 06 

-0.3392 xlO" 06 

-0.1190xl0 -06 

0.1093 xlO’ 06 

10 

8 

0.1427 xlO' 06 

-0.2123 xlO" 06 

0.4874xl0-° 7 

-0.7202 xlO" 07 

12 

8 

0.3999 xlO" 07 

-0.6142 xlO" 07 



14 

8 

0.6291 xlO" 08 

-0.6646 xlO" 07 


1 

16 

8 

-0.1185xl0- 07 

-0.3751 xlO" 07 

-0.3917xl0-° 7 

0.3830 xlO" 08 

18 

8 

0.5125 xlO" 07 

-0.2834 xlO’ 07 



20 

8 

0.1794xl0-° 7 

-0.1104xl0-° 7 

0.4059 xlO" 08 

0.1485 XlO- 07 

22 

8 

-0.2004xl0 - ° 7 

-0.1722xl0 _o7 

-0.3188xl0-° 7 

0.2914xl0-° 8 

24 

8 

0.1932 xlO" 07 

-0.2081 xlO" 07 



26 

8 

0.7610xl0-° 8 

-0.1681 xlO" 07 

0.9028 xlO" 09 

-0.3295 XlO- 08 


8 

0.6911 xlO" 08 

-0.1509 xlO" 07 

g‘ffTTTTT*toi 


30 

8 

0.3667 xlO -08 

0.2371 xlO" 09 



32 


0.1242 xlO" 07 

-0.2447 xlO" 08 

0.9457 xlO" 08 

0.7227xl0 -08 

34 

8 

-0.7679 xlO" 08 

-0.5526 xlO’ 08 

B‘ffTlTITtVm 


36 

8 

0.3052 xlO" 08 

-0.1028 xlO" 07 

0.9438 xlO' 10 

-0.3564X10" 08 

38 

H 

0.5463 xlO" 08 

-0.5407 xlO" 08 



40 

H 

0.7907xl0-° 8 

-0.2545 xlO" 08 



42 

8 

0.4186 xlO" 08 

-0.8110 xlO" 08 



44 

8 

-0.6108 xlO" 08 

-0.5257 xlO -08 

-0.8194xl0-° 8 

-0.1681 xlO" 09 

46 

8 

0.5875 xlO" 09 

-0.1010 xlO" 08 

-0.3281 xlO" 09 

0.3699 xlO" 08 

48 

8 

0.3773xl0 -08 

-0.2395 xlO" 08 

0.2099 xlO" 08 

0.2050 xlO" 08 

50 

8 

-0.3182xl0-° 8 

-0.6258 xlO" 08 

-0.4874xl0-° 8 

-0.2563 xlO" 08 

52 

8 

0.1505xl0-° 8 

-0.3795xl0-° 8 

0.2233 xlO" 09 

-0.3470 xlO” 09 

54 

8 

0.2067 xlO" 08 

-0.5294 xlO" 08 

0.9532 xlO" 09 

-0.1832 xlO' 08 

56 

8 

-0.7950 xlO" 09 

-0.9549 xlO' 08 

-0.1446 xlO" 08 

-0.6491 xlO“ 08 

58 

8 

0.2954xl0-° 8 

0.1545 xlO' 08 

0.2019 xlO' 08 

0.4413 xlO" 08 

60 

8 

0.4263 xlO" 08 

-0.2996 xlO" 08 

0.3292 xlO" 08 

0.4599 xlO- 10 


Table 18: Recovered Coefficients for Block m — 8 , [ / (2 = 0 
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estimate 

truth 

l 

m 

c lm 

Sim 

Clm 

Sim 

IE9 


-0.1531 xlO" 06 

-0.3231 xlO" 07 

0.1243 xlO" 06 

-0.4088 xlO" 07 

12 

9 

-0.1415 xlO" 06 

0.1643 xlO' 07 

0.3980 xlO" 07 

0.3063 xlO" 07 

14 

9 

-0.1058 xlO' 06 

0.1591 xlO" 07 

0.1821 xlO -07 

0.2687 xlO" 07 

16 

9 

-0.1099 xlO" 06 

-0.5081 xlO" 07 

-0.2389 xlO- 07 

-0.5027 xlO" 07 

18 

9 

-0.8810 xlO" 07 

0.1652 xlO" 07 

-0.2129 xlO" 07 

0.2512 xlO" 07 

20 

9 

gfewsiiiaiifla 


0.1320 xlO' 07 


22 

9 



0.4872 xlO" 08 

9 



EEIllQl^B 


-0.1182 xlO" 07 

-0.2214X10" 07 



-0.4289 xlO“ 07 

0.1558xl0-° 8 

-0.9284xl0-° 8 

0.4153 xlO" 08 





0.8193 xlO" 08 

-0.4057 xlO” 08 

30 

9 

-0.3196 xlO" 07 

-0.5596 xlO" 08 

-0.7258 xlO~ 08 

-0.3901 xlO" 08 

32 

9 

-0.2041 xlO" 07 

-0.3481 xlO" 08 

0.2354 xlO" 08 

-0.1813xl0-° 8 

34 

9 

-0.2003 xlO" 07 

-0.8035 xlO' 09 

0.4458 xlO" 09 

0.8343 xlO- 09 

36 

9 

-0.1757 xlO" 07 

-0.1722 xlO" 08 

0.1012 xlO' 08 

-0.1510 xlO" 09 

38 

9 

-0.1060 xlO' 07 

-0.3810 xlO' 08 

0.6977 xlO" 08 

-0.1801 xlO’ 08 

40 

9 

-0.1492 xlO" 07 

0.1750 xlO' 09 

-0.1622 xlO' 09 

0.1263 xlO" 08 

42 

9 

-0.1648 xlO" 07 

0.1060 xlO" 08 

-0.1354xl0-° 8 

0.1937xl0-° 8 

44 

9 

-0.1347 xlO" 07 

-0.9001 xlO" 08 

0.2005 xlO" 09 

-0.7535 xlO" 08 

46 

9 

-0.7349 xlO" 08 

-0.5857 xlO" 09 

0.5190 xlO" 08 

0.1447 xlO" 08 

48 

9 

-0.1586 xlO" 07 

0.2641 xlO" 08 

-0.3948 xlO" 08 

0.3988 xlO" 08 

50 

9 

-0.1370 xlO' 07 

0.4798 xlO" 09 

-0.2317 xlO" 08 

0.9108 xlO- 09 

52 

9 

-0.1316 xlO* 07 

-0.4716 xlO" 08 

-0.2667 xlO~ 08 

-0.3721 xlO' 08 

54 

9 

-0.8838 xlO" 08 

0.4068 xlO" 08 

0.1330 xlO" 08 

0.6034X10" 08 

56 

9 

-0.7253 xlO” 08 

0.2354 xlO" 08 

0.2965 xlO" 08 

0.3267 xlO' 08 

58 

9 

-0.1243 xlO" 07 

-0.4833 xlO" 08 

-0.3562 xlO" 08 

-0.3847 xlO' 08 

60 

9 

-0.7154xl0- 08 

0.7376 xlO' 09 

0.2488 xlO" 08 

0.2052 xlO" 08 


Table 20: Recovered Coefficients for Block m = 9 , |/|2 = 0 
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estimate 

truth 

l 

TO 

c lm 

Si m 

Clm 

Sim 

9 

9 

0.2161 xlO- 06 

0.1 167 xlO" 06 

-0.4490 xlO" 07 

0.8300 xlO" 07 

11 

9 

0.9118 xlO- 07 

0.4504 xlO" 07 

-0.3117xl0-° 7 

0.5650 xlO- 07 

13 

9 

0.6566 xlO” 07 

0.3984xl0-° 7 


hT^TTHTVMI 

15 

9 

0.5580 xlO” 07 

0.2834xl0-° 7 



17 

9 

0.3260 xlO' 07 

-0.3917xl0-° 7 



19 

9 

0.3127 xlO" 07 

-0.5876 xlO" 08 

0.7484x10-°® 

-0.3321x10-" 

21 

9 

0.3447 xlO" 07 

0.1064xl0-° 7 

0.1455xl0-° 7 

0.1543 xlO” 07 

23 

9 

0.1844xl0-° 7 

-0.1380 xlO" 07 

0.8397xl0-° 9 

-0.1206 xlO" 07 

25 

9 

-0.1195 xlO" 07 

0.6443 xlO” 08 

-0.2786 xlO- 07 

0.1184xl0-° 7 

27 

9 


BitnwBii 


HfTlnfTTFWi 

29 

9 





31 

9 

0.7488 xlO" 08 

-0.5687 xlO' 09 

-0.1704xl0- 08 

0.2351x10-" 

33 

9 

0.7192 xlO" 08 

-0.5837 xlO" 09 

-0.1606 xlO" 08 

0.1999x10-" 

35 

9 

0.7294xl0-° 8 

-0.2835 xlO" 08 

-0.3823x10-" 

-0.8356x10-" 

37 

9 

0.7441 xlO" 08 

-0.6238 xlO’ 08 

0.5915x10-" 

-0.4685x10-" 

39 

9 

0.1418 xlO” 07 

0.4807 xlO- 08 

0.7844x10-" 

0.7154x10-" 

41 

9 

-0.3601 xlO- 09 

0.2978 xlO' 08 

-0.6638x10-" 

0.4913x10-" 

43 

9 

0.5809 xlO" 08 

-0.8692 xlO" 08 

0.6125x10-" 

-0.8150x10-" 

45 

9 

0.1007 xlO" 07 

-0.6726 xlO' 08 

0.5668x10-" 

-0.4822x10-" 

47 

9 

0.3694xl0 -08 

0.2529 xlO" 08 

-0.1809x10-" 

0.4123x10-" 

49 

9 

0.3780 xlO" 08 

0.5580 xlO" 08 

-0.2210x10-" 

0.6961x10-" 

51 

9 

0.5860 xlO" 08 

-0.3246 xlO" 08 

0.2005x10-" 

-0.2567x10-" 

53 

9 

0.6105 xlO" 08 

-0.3452 xlO" 08 

0.2767x10-" 

-0.1838x10-" 

55 

9 

0.4237 xlO" 08 

0.2124xl0 -08 

0.4647x10-" 

0.2930x10-" 

57 

9 

0.5303 xlO" 08 

-0.2800 xlO" 08 

0.1317x10-" 

-0.2636x10-" 

59 

9 

0.2607 xlO" 08 

-0.6082 xlO" 09 

-0.3646x10-" 

0.1065x10-" 


Table 21: Recovered Coefficients for Block m = 9 , |/| 2 = 1 
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estimate 

truth 

l 

m 

c, m 


C im 

■! Sim 

10 

10 

-0.5140 xlO" 06 

-0.4307xl0 -07 

0.1075 xlO" 06 

-0.2786 xlO" 07 

12 

10 

-0.2473 xlO" 06 

-0.1140xl0- 07 

bhbh 

EHH&BS9 

14 

10 

-0.1262 xlO" 06 

-0.3844xl0 -07 



16 

10 

-0.8863 xlO" 07 

-0.2367 xlO" 07 

0.1315X10" 08 

0.5816 xlO" 08 

18 

10 

-0.3985 xlO" 07 

-0.2389 xlO" 07 



20 

10 

-0.6424 xlO" 07 

-0.2814xl0-° 7 

-0.2261 xlO" 07 

-0.6292 xlO" 08 

22 

10 

-0.3040 xlO’ 07 

0.1949 xlO" 08 

0.2076 xlO" 08 

0.2270 xlO" 07 

24 

10 

-0.1668 xlO" 07 

0.6901 xlO" 08 

0.7879 xlO- 08 

0.2529 xlO" 07 

26 




-0.1040 xlO" 07 

-0.1149xl0-° 7 





-0.8553 xlO" 08 

0.1151xl0-° 7 

30 

10 

-0.1347 xlO" 07 

-0.1563 xlO' 07 

-0.1655 xlO’ 08 

-0.4794xl0- 08 

32 

10 

-0.9087 xlO" 08 

-0.1613 xlO" 07 

0.1360 xlO" 08 

-0.6172 xlO" 08 

34 

10 

-0.1675 xlO" 07 

-0.1073 xlO" 07 

-0.8343 xlO" 08 

-0.1238 xlO- 08 

36 

10 

-0.8031 xlO" 08 

-0.3227xl0-° 8 

-0.2170 xlO" 09 

0.5186xl0-° 8 

38 

10 

-0.9255 xlO" 08 

-0.1178xl0-° 7 

-0.3025 xlO" 08 

-0.4598 xlO" 08 

40 

10 

-0.1060 xlO" 07 

0.1348 xlO" 08 

-0.4826 xlO" 08 

0.8558 xlO- 08 

42 

10 

-0.9558 xlO" 09 

-0.3201 xlO" 08 

0.4033 xlO" 08 

0.3632 xlO" 08 

44 

10 

-0.8301 xlO" 08 

-0.9874 xlO- 08 

-0.3711 xlO' 08 

-0.4245 xlO" 08 

46 

10 

-0.6224xl0- 09 

-0.5392 xlO" 08 

0.3573xl0-° 8 

0.4634xl0-° 9 

48 

10 

-0.3753 xlO" 08 

-0.1804xl0-° 8 

-0.4231 xlO" 09 

0.3661 xlO" 08 

50 

10 

-0.7997 xlO" 08 

-0.4131 xlO’ 08 

-0.4087xl0- 08 

0.7728 xlO' 09 

52 

10 

-0.3065 xlO" 08 

-0.5195 xlO" 08 

0.1736 xlO' 09 

0.7139xlO' 10 

54 

10 

-0.2656 xlO" 09 

-0.3731 xlO" 08 

0.1743xlO -08 

0.1307xl0-° 8 

56 

10 

-0.5989 xlO" 08 

-0.4236 xlO" 08 

-0.3497 xlO” 08 

-0.2255 xlO" 09 

58 

10 

-0.3760 xlO" 08 

-0.7944 xlO" 08 

-0.8522 xlO" 09 

-0.4618xl0-° 8 

60 

10 

-0.3664xl0-° 8 

-0.5742 xlO" 08 

-0.1824 xlO" 08 

-0.1196xl0- 08 


Table 22: Recovered Coefficients for Block m = 10 , |/j 2 = 0 
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estimate 

truth 



c, m 

Sim 

Clm 

S',™ 

11 

10 

0.3627 xlO -07 

0.5296 xlO' 07 

-0.4870 xlO" 07 

-0.1692 xlO" 07 

13 

10 

0.8194xl0~ O7 

0.2467 xlO" 07 

0.3043 xlO" 07 

-0.4290 xlO" 07 

E3 

m 

0.5272 xlO" 07 

0.5619 xlO" 07 

0.5619 xlO” 08 

0.6209 xlO" 08 

EH 

E9 

0.4956 xlO- 07 

0.5642 xlO~ 07 

0.7784xl0-° 8 

0.1316xl0-° 7 

19 

10 

0.3049 xlO" 07 

0.3325 xlO" 07 

-0.1068 xlO" 07 

-0.7629 xlO" 08 

21 

10 

0.2775 xlO" 07 

0.3662 xlO" 07 

-0.7907 xlO" 08 

0.6094xl0- xo 



0.4146 xlO- 07 

0.3367 xlO -07 

0.8725 xlO” 08 

-0.1901 xlO" 08 

El 


0.3967 xlO -07 

0.2854 xlO" 07 

0.8381 xlO" 08 

-0.6219 xlO" 08 



0.1951 xlO’ 07 

0.3215 xlO" 07 

-0.1019xl0- 07 

0.2050 xlO" 08 

29 

10 

0.3267xl0-° 7 

0.3211 xlO- 07 

0.5527xl0-° 8 

0.1588 xlO" 08 

31 

10 

0.2386 xlO” 07 

0.2500 xlO' 07 

-0.2927 xlO" 08 

-0.3669 xlO" 08 

33 

10 

0.2167xl0-° 7 

0.2401 xlO" 07 

-0.4557 xlO" 08 

-0.3459 xlO" 08 

35 

10 

0.2075 xlO” 07 

0.3120 xlO" 07 

-0.3516 xlO" 08 

0.4851 xlO" 08 

37 

10 

0.2389 xlO” 07 

0.2981 xlO" 07 

-0.3048 xlO" 10 

0.3140 xlO" 08 

39 

10 

0.2051 xlO” 07 

0.2972 xlO" 07 

-0.2746 xlO' 08 

0.4983 xlO" 08 

41 

10 



0.3415X10- 08 

-0.2164xl0-° 9 

43 

10 

0.1964xl0-° 7 

0.2702 xlO" 07 

-0.3208 xlO" 08 

0.1770 xlO" 08 

45 

10 

0.2347 xlO" 07 

0.2397 xlO" 07 

0.2470 xlO- 08 

-0.2523 xlO" 09 

47 

10 

0.2321 xlO’ 07 

0.2615 xlO" 07 

0.1688 xlO" 08 

0.1998 xlO" 08 

49 

10 

0.1524xlO -07 

0.2509 xlO" 07 

-0.6413 xlO" 08 

0.6391 xlO" 09 

51 

10 

0.2268 xlO' 07 

0.2000 xlO" 07 

0.2505 xlO- 08 

-0.4520 xlO" 08 

53 

10 

0.2870 xlO’ 07 

0.2261 xlO" 07 

0.8869 xlO" 08 

-0.1494xl0-° 8 

55 

10 

0.1950xl0-° 7 

0.2777 xlO" 07 

-0.1480 xlO" 08 

0.4015 xlO" 08 

57 

10 

0.1788 xlO" 07 

0.2602 xlO” 07 

-0.2301 xlO" 08 

0.2124X10- 08 

59 

10 

0.2292 xlO- 07 

0.1963 xlO" 07 

0.2990X10- 08 

-0.4309 xlO" 08 


Table 23: Recovered Coefficients for Block m = 10 , |/|2 = 1 
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Figure 5: Observed SST Range Rate Perturbation Spectral Amplitude: m — 2, |/| 2 = 0 
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Figure 6: Observed SST Range Rate Perturbation Spectral Phase: m = 2, |Z| 2 = 0 


60 



-100000. -80000. -60000. -40000. -20000. .. 20000. 40000. 60000. 80000. 100000. 


Figure 7: Observed SST Range Rate Perturbation Spectral Amplitude: m — 2, |/|2 = 1 
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Figure 8: Observed SST Range Rate Perturbation Spectral Amplitude: m — 3, |/|2 = 0 
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Figure 9: Observed SST Range Rate Perturbation Spectral Amplitude: m = 3, |Z |2 = 1 
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Spectral Amplitude: m = 4, |/| 2 = 0 
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Figure 11: Observed SST Range Rate Perturbation Spectral Amplitude: m = 4, |/| 2 = 1 
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Figure 12: Observed SST Range Rate Perturbation Spectral Amplitude: m = 5, | / j 2 = 0 
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Figure 13: Observed SST Range Rate Perturbation Spectral Amplitude: m = 5, j/j 2 = 1 
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Figure 14: Observed SST Range Rate Perturbation Spectral Amplitude: m = 6, |/| 2 = 0 
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Figure 15: Observed SST Range Rate Perturbation Spectral Amplitude: m = 6, |/|2 = 1 
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Figure 19: Observed SST Range Rate Perturbation Spectral Amplitude: m = 8, |/|2 = 1 
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Figure 20: Observed SST Range Rate Perturbation Spectral Amplitude: m = 9, |/| 2 = 0 
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Figure 25: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 1, |/| 2 = 1 
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Figure 27: Predicted SST Range Rate Perturbation Spectral Phase: m = 2, |/j 2 = 0 
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Figure 28: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 2, j 2 = 1 
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Figure 29: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 3, | / J 2 = 0 
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Figure 30: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 3, |/| 2 = 1 
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Spectral Amplitude: m — 5, |/| 2 = 1 
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Figure 35: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 6, |/|2 = 0 
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Figure 37: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 7, \l\i = 0 
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Figure 38: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 7, \l \2 = 1 
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Figure 39: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 8, |/| 2 = 0 
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Figure 41: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 9, | Z I 2 = 0 
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Figure 42: Predicted SST Range Rate Perturbation Spectral Amplitude: m = 9, |/| 2 = 1 
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Figure 44: Predicted SST Range Rate Perturbation Spectral Amplitude: m — 10, |Z |2 = 1 
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